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ABSTRACT 

A  two  dimensional,  incompressible  viscous  inviscid  interaction  computer  code,  de- 
signed to  compute  cascade  flows,  was  investigated.  Comparison  of  the  flow  character- 
istics predicted  by  the  code  with  experimentally  available  data  indicates  that  the  code 
predicts  reasonably  well  flow  parameters  on  lightly  loaded  cascades.  However,  the  code 
fails  to  predict  correctly  the  actual  boundary  layer  development  and  the  velocity  dis- 
tribution for  highly  loaded  cascades.  It  is  concluded  that  further  improvement  of  the 
code  is  needed  and  recommendations  are  presented  to  achieve  the  required  improve- 
ments. 
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I.     INTRODUCTION 

The  need  for  better  and  more  efficient  gas  turbines  requires  the  availability  of  cheap 
and  reliable  design  tools  for  blades  used  in  compressors  and  turbines.  Computational 
methods  are  the  preferred  choice  for  such  a  design  tool,  considering  the  cost  and  com- 
plexity of  wind  tunnel  experiments. 

Among  the  computational  methods  available  today,  the  logical  choice  seems  to  be 
a  computer  code  that  can  solve  directly  the  full  Navier  Stokes  equations.  However,  given 
the  state  of  the  art  in  both  algorithms  and  computer  hardware,  such  Navier  Stokes 
solvers  are  restricted  only  to  supercomputers,  and  even  then  the  computation  time  is 
quite  long. 

In  order  to  enable  fast  and  efficient  computations,  the  viscous  inviscid  interaction 
code  was  developed  by  Cebeci  [Ref.  1].  The  approach  used  in  this  code  is  to  solve  the 
outer  flow  field  using  potential  methods,  and  solving  the  boundary  layer  flow  using  a 
boundary  layer  method  subject  to  an  interaction  law,  that  couples  the  inner  and  the 
outer  flows.  This  interaction  law  is  needed  because  classical  boundary  layer  methods  fail 
in  areas  of  flow  reversal  and  separation,  which  are  very  common  in  real  life  flows. 

The  viscous  inviscid  interaction  code  was  originally  developed,  and  successfully  used, 
for  flows  about  single  airfoils.  It  was  later  adapted  to  cascade  flows. 

In  this  thesis  the  applicability  of  the  code  to  cascades  was  investigated  by  comparing 
its  results  to  experimental  data.  It  was  found  that  although  the  code  can  reasonably 
predict  experimental  results  in  some  cases,  it  still  needs  improvements  before  it  can  be 
applied  generally  as  a  reliable  design  tool. 

A  major  restriction  in  improving  the  code  is  the  lack  of  a  wide  data  base  of  appro- 
priate experimental  results.  Some  of  the  key  elements  in  the  code,  like  transition  and 
turbulence  modelling,  are  based  on  empirical  correlations,  and  more  detailed  and  accu- 
rate experiments  should  be  performed,  before  a  better  understanding  of  these  phenom- 
ena can  be  achieved. 

In  the  following,  the  theoretical  background  of  the  code  is  presented  in  Chapter  II. 
a  description  of  the  code  in  Chapter  III,  the  results  and  discussions  are  presented  in 
Chapter  IV  and  the  conclusions  and  recommendations  in  Chapter  V.  A  listing  of  the 
computer  program  is  given  in  Appendix  A. 


II.     CASCADE  FLOW  PROBLEM  FORMULATION 

This  chapter  outlines  the  theoretical  background  of  the  viscous  inviscid  interactive 
method  used  in  the  computer  code  to  investigate  cascade  flows.  Only  the  major  steps  in 
the  mathematical  developments  will  be  described  here.  A  detailed  description  of  the 
theory  and  the  numerical  methods  is  given  by    Cebeci  and  Bradshaw  [Ref.  2]  and  by 
Krainer  [Ref.  3]  on  which  this  chapter  is  based. 

A.     INVISCID  FLOW  METHOD 

Inviscid  flow  is  the  first  building  block  of  the  flow  and  is  solved  using  the  panel 
method.  The  incompressible  two  dimensional  outer  flow  must  satisfy  the  Laplace 
equation: 

V2cD  =  0  , 

subject  to  the  boundary  conditions  on  the  surface  of  the  blade: 


en 


where  the  commonly  used  boundary  condition  of  zero  normal  velocity  on  the  surface  is 
replaced  by  a  specified  blowing  velocity  vw  to  allow  for  the  effect  of  the  boundary  layer 
on  the  outer  flow. 

In  addition,  the  Kutta  condition  must  be  satisfied,  in  order  to  prevent  the  existence 
of  discontinuous  pressure  distribution  near  the  trailing  edge  of  the  blade. 

Since  the  Laplace  equation  is  linear,  a  solution  to  a  complex  flow  field  can  be  built 
by  superposition  of  solutions  of  elementary  flows.  The  elementary  flows  used  in  the 
panel  method  are  the  uniform  parallel  flow  and  flows  about  singularities  (sources  and 
vortices). 

The  panel  method  is  based  on  replacing  each  blade  by  a  distribution  of  sources  and 
vortices  on  its  surface.  The  surface  is  divided  into  a  finite  number  of  straight  segments, 
called  panels. 

On  each  panel,  a  uniform  source  distribution  and  a  uniform  vorticity  distribution  is 
assumed.  The  source  strength  at  each  panel  is  set  to  satisfy  the  boundary  condition  at 
the  midpoint  of  the  panel  (called  the  control  point),  and  so,  in  general  the  source 


strength  will  vary  from  panel  to  panel.  The  vorticity  strength  is  assumed  to  be  the  same 
for  all  the  panels  and  is  set  to  satisfy  the  Kutta  condition. 

The  cascade  is  defined  as  an  infinite  row  of  similar  blades,  each  one  modelled  by 
panels  of  source  and  vortex  distributions.  The  flow  at  each  point  is  found  by  summing 
the  contributions  of  all  the  singularities  on  all  the  blades  ,  and  the  uniform  parallel  flow. 

A  useful  concept  in  dealing  with  such  flows  is  the  concept  of  influence  coefficients. 
An  influence  coefficient  is  defined  as  the  velocity  at  a  point  induced  by  a  unit  strength 
singularity  placed  at  some  other  field  point.  The  influence  coefficients  are  a  function  of 
geometry  and  so  can  be  computed  for  a  given  cascade  and  a  given  choice  of  panel  ge- 
ometry. 

Usin2  the  influence  coefficients,  the  normal  and  tangential  velocities  at  each  control 
point  can  be  written  as  a  function  of  the  unknown  source  strength  of  each  panel  and  the 
unknown  vortex  strength.  Using  the  boundary  conditions,  by  equating  the  normal  ve- 
locity at  each  control  point  to  the  prescribed  blowing  velocity  v„,  and  using  the  Kutta 
condition  (which  requires  equal  velocities  on  the  upper  and  on  the  lower  panels  at  the 
trailing  edge),  a  system  of  linear  equations  is  constructed. 

By  solving  the  system  of  equations,  the  strength  of  the  sources  and  vortices  is  found, 
and  the  velocities  (and  the  pressures)  can  be  computed  everywhere  in  the  flow  field. 

The  velocity  distribution  on  the  surface  of  the  blade,  computed  by  the  panel  method, 
is  used  as  the  boundary  condition  for  the  boundary  layer  flow  calculations. 

It  should  be  noted  that  in  the  panel  method  as  used  in  the  present  computer  code, 
there  is  no  modelling  of  the  wake,  and  its  effect  on  the  flow  field  is  ignored. 

B.     VISCOUS  FLOW  METHOD 

Viscous  flow  is  the  second  building  block  of  the  cascade  flow  and  it  is  applied  to  the 
thin  boundary  layer  near  the  blade  surface. 
1.     Boundary  Layer  Theory 

The  boundary  layer  concept,  first  suggested  by  L.  Prandtl,  assumes  that  the  flow 
field  can  be  divided  into  an  outer  flow  where  the  viscous  effects  are  negligible  compared 
to  inertia  effects,  and  a  thin  layer  close  to  the  surface  where  the  viscous  effects  cannot 
be  neglected.  The  complete  presentation  of  the  boundary  layer  theory,  and  the  devel- 
opment of  the  boundary  layer  equations,  is  given  by  Schlichting  [Ref.  4]. 

Under  the  assumptions  of  two  dimensional  incompressible  thin  boundary  layer 
flow,  the  Navier  Stokes  equations  and  the  continuity  equation  reduce  to: 


-SSL  +  -2L   =  o 
ex        cy 


w  -r^1    h  v  -~  =  w.  — — —  v 


5jc  cy        e  <£c  6y  \     8y  J  ' 

with  the  boundary  conditions: 

v  =  0  u(x,0)  =  0  ,  v(x,0)  =  0  , 

y=ye  "(X&'e)  =  ue{x), 

V, 

where  b  denotes  1  +  — . 

Writing  the  velocity  components  in  terms  of  a  stream  function  *F 
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This  eliminates  the  continuity  equation  (which  the  stream  function  satisfies  by  defi- 
nition). 

Introducing  the  Falkner  Skan  transformation: 


j{xy)  =       . *l*{xy), 

yjUeVX 


the  momentum  equation  and  the  boundary  conditions  transform  to: 

(¥")'  +  2n4±ff"  +  m[l  -  (/")2]  =  x(f  4~  -/"  -|f  )  , 

^  \  OX  OX     J 


V  =  0       f'(x,0)  =  0  ,J{xy)  =  o , 
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n  =  ne  fix,  rje)  =  l , 

where  m  is  defined  by: 


m  = 


x    due 


ue    dx 


The  third  order  differential  equation  can  be  reduced  to  a  system  of  first  order 
differential  equations  by  introduction  of  two  new  variables  U  and  V  : 

U  -  /'. 

V  -  L'  , 

(bY)'+-2^r-JV  +  m{\-U2)  =  x(\J-^--y4-), 
2  \      ox   ■        ox  J 

with  the  boundary  conditions: 

n  =  o       u(.v,o)  =  o  ,j[x,o)  =  o , 

The  next  step  is  to  use  a  finite  difference  approach  to  solve  the  equations.  The 
box  method  is  applied  using  central  differencing  in  both  the  x  and  rj  directions,  and 
satisfying  the  equations  midway  between  nodes. 

Applying  the  box  method  results  in  a  system  of  nonlinear  equations  in  the  un- 
known variables  (which  are  f,  U  and  V  in  each  node  along  the  rj  direction  at  the  current 
x  station). 

In  order  to  solve  the  nonlinear  system  the  Newton  iterative  procedure  is  used, 
linearizing  the  equations  first  about  the  solution  at  the  adjacent  upstream  station,  and 
then  about  the  preceding  iteration.   The  linearization  is  performed  by  letting: 

fj  =fj  +%■  , 

t  t/,  K  t  •(',  K—  1       ,        rT  •/,  K 

Lj      =    Lj         +  0  L.J     , 


\  rl,  K  »  A,  K—  1       .        rir/.IC 


where: 

•  /  denotes  location  in  the  x  direction 

•  j  denotes  location  in  the  y  (  r\  )  direction 

•  k  indicates  the  iteration  counter 

This  linearization  results  in  a  system  of  linear  equations  for  the  unknown  in- 
crements: df/x ,  SL';K  and  <5V;-r . 

This  system  of  equations  is  solved  repeatedly  until  the  changes  in  the  unknowns 
are  small  enough.  Since  the  system  is  block  tridiagonal,  Keller's  block  elimination 
method  is  used. 

The  method  described  so  far,  is  a  direct  boundary  layer  method.  It  can  be  used 
as  long  as  the  flow  does  not  separate.  Whenever  separation  or  flow  reversal  occurs,  and 
a  zero  skin  friction  coefficient  is  encountered,  the  equations  become  singular  and  the 
calculations  will  break  down. 

2.     Interactive  Boundary  Layer  Method 

The  interactive  boundary  layer  method  is  designed  to  overcome  the  difficulties 
encountered  at  regions  of  flow  reversal  and  separations.  In  such  areas  the  external  ve- 
locity is  substantially  changed  by  the  viscous  effects  and  can  no  longer  be  considered  as 
a  known  boundary  condition  for  the  boundary  layer  flow. 

The  general  approach  to  the  solution  is  the  same  as  for  the  direct  method  but. 
since  the  outer  flow  is  unknown,  the  velocity  at  the  edge  of  the  boundary  layer  is  written 
as: 


u(x,ye)  =  uel  +  — 


where: 

1.  ue(x,ye)    is  the  total  velocity  at  the  edge  of  the  boundary  layer. 

2.  uel(x)    is  the  velocity  as  computed  by  the  inviscid  method. 
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— -  (ueS') is  the  Hilbert  integral. 
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The  numerical  solution  of  the  boundary  layer  equations  follows  the  same  steps 
as  for  the  direct  method,  but  with  some  changes. 

The  transformation  of  the  stream  function  and  they  coordinate  uses  a  constant 
velocity  u0  as  a  scaling  factor,  and  a  scaled  velocity  w  is  introduced: 


/  "o 


fx,  rj)  =        . \j/(x,y)  , 

JUqVX 


ue{x.y) 
w  = 


u. 


0 


Using  this  transformation,  the  boundary  layer  equations  become  a  system  of 
first  order  differential  equations: 

/'  =  U , 
U'  =  V, 

W  =  0, 

(bvy  +  ±-jv  +  x\v^-  =  x(v-^-y^), 

2  ex  \       ex  ex  J 

with  the  boundary  conditions: 

n  =  o       u(jr,0)  -  o  ,y(jcfo)  =  o , 

q  =  r\e  U{x,  t]e)  =  W(jf,  vie)  • 


w(*.  *.)  =    ^  +  T  J  ^  (\/t£  W«.  *X  -JK,  ^     ^ 
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The  finite  difference  box  method  is  used  to  solve  the  equations,  in  the  same  way 
as  it  was  used  for  the  direct  case,  but  with   two  additions: 


1.  In  areas  of  flow  reversal  the  term  ucu'cx  is  omitted  to  assure  stable  integration 
(the  FLARE  approximation). 

2.  The  edge  velocity,  Wj  (where  J  denotes  the  edge  station)  which  involves  integration, 
is  approximated  by  : 

W,  -  g{  +  CaCWjjf,  -A)  , 

where  g,  and  c„  are  obtained  from  the  numerical  approximation  to  the  Hilbert  in- 
tegral (which  will  be  presented  in  the  next  section  ). 

By  using  central  diffrencing  to  approximate  the  differential  equations,  a  system 
of  nonlinear  algebraic  equations  is  obtained  for  the  unknown  variables  (which  are 
fj ,  Uj ,  V;  and  Wj).  To  solve  the  system  of  equations,  the  system  is  linearized  by  the 
Newton  iterative  procedure,  and  the  resulting  linear  system  is  solved  (for  the  new  un- 
known variables  which  are  the  increments  8f>* ,  dV/K ,  SV'/K  and  5W'f" ). 

The  solution  of  the  system  is  repeated  until  the  change  in  the  increments  is 
negligible  compared  to  the  preceding  iteration,  and  the  whole  process  is  performed  again 
at  the  next  downstream  station. 
3.     Interactive  Model 

The  interactive  model  is  used  to  couple  the  boundary  layer  to  the  external  flow. 
It  is  needed  in  areas  where  strong  interaction  occurs,  and  both  the  boundary  layer  and 
the  outer  flow  must  be  solved  simultaneously.  The  interaction  model  provides  the  outer 
boundary  condition  to  the  boundary  layer  calculations  by  adding  a  correction  term  to 
the  external  velocity  computed  by  the  inviscid  flow  method. 

The  external  velocity  is  assumed  to  consist  of  a  potential  flow  term  (  uei{x)  )  and 
a  correction  term  due  to  viscous  effects  (  uei  (x)  ): 

"eW  =   uel{x)  +   ue6{x). 

The  viscous  effect  is  obtained  by  a  surface  distribution  of  sources  on  the  blade  (a  con- 
cept first  suggested  by  Lighthill  [Ref.  5]).  The  normal  velocities  at  the  surface  of  the 
blade,  induced  by  these  sources,  displace  the  streamlines  from  the  surface  in  the  same 
way  that  the  actual  boundary  layer  displaces  them: 

dS'(x)  v(x,  d*) 


dx  ue(x)     ' 

Where  \(x,  S')  is  the  normal  velocity  at  the  displaced  surface. 


Assuming  that  the  surface  can  be  approximated  to  he  a  flat  plate,  the  normal 
velocity  will  be  half  the  local  source  strength  o  (x).  Assuming  also  that  the  inviscid  ve- 
locity doe^  not  chance  across  the  boundarv  laver.  the  local  source  strength  will  be: 


a(x) 


=  v{x,0)  =  v{x,  S  ) 


"0 


cv 


cy 


dy 


d 

dx 


iVe*  ) 


The  local  horizontal  velocity  induced  by  the  source  distribution,  is  the  cor- 
rection term  to  the  inviscid  velocity,  and  can  be  represented  by  the  Hilbert  integral: 


rxh 


r\l)    r_   =    J_ 


d 
dl 


[u0b 


dl 


x 


The  integration  is  carried  out  on  all  the  sources  on  the  surface,  since  the  horizontal 
velocity  is  influenced  by  all  the  sources. 

The  Hilbert  integral  is  then  approximated  by  a  finite  series: 


71 


d    ,     ~\     dl 


dl 


x  —  c 


cik{ue6  )  - . 


/:=! 


Where  cih  is  a  matrix  of  interaction  coeflicients  which  are  functions  of  the  geometry  only 
(  i  denotes  the  chordwise  position  where  uei  is  evaluated  and  k  is  the  location  of  the 
source  which  effects  ue6  ). 

Since  the  computation  of  uei  involves  values  of  3'  downstream  of  the  current  x 
location,  which  are  not  known  yet,  these  terms  are  taken  from  the  previous  iteration 
using  a  relaxation  formula. 
4.     Turbulence  Model 

The  turbulence  model  used  here  is  the  algebraic  eddy  viscosity  formulation  of 
Cebeci  and  Smith  [Ref.  6].  According  to  the  model  used  in  the  present  computer  code. 
the  eddy  viscosity  v,  is  defined  by  two  different  expressions,  for  the  inner  region  and  for 
the  outer  region: 


v,=  Wl-exp(-iy 


8y 


Yxr 


for    0<y<yc, 


/A 


for  yc  <y  <  S 


Where: 


A  = 


26v 

C?/     \l/2 


1  +  5.5(W<))6 


a  = 


0.0168 


1  -  p 


cu.'cv 


2.5 


P    = 


1  +  2RT(2  -  RT) 


for     RT  <  1  , 


/?  = 


1  +  RT 

Rt 


for      RT  >  1  , 


RT  = 


(  -  u'v') 


max 


The  distance  from  the  wall  to  the  point  between  the  two  regions,  ye,  is  chosen  such  that 
the  viscosity  will  be  continuous. 
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The  intermittency  factor,  yf,   is  defined  by: 


ytr  =  1  -  exp 


3 

Gyv 


rx 


di 


Where: 

•  R,r,  is  the  Reynolds  number  based  on  external  velocity  and  transition  location. 

•  Gv  is  an  empirical  constant,  originally  assigned  the  value  1200. 


Cebeci  and  Bradshaw  [Ref.  2,  p. 246]  described  a  different  expression  for  the 
variable  A  in  the  inner  region  viscosity  formula: 


A  = 


26* 


(i  -  n.zp-)]!2(v^}12 

\      °y   /max 


Where: 


via, 


P      - 


v_£«.\3/2 
.    oy 


dx 


This  version  of  the  turbulence  model  was  not  implemented  in  the  original  computer 
code.  During  the  work  on  this  thesis,  the  effect  of  the  modified  turbulence  model  was 
investigated. 

A  different  intermittency  distribution  was  implemented  successfully  by  Rodi  and 
Schonung  [Ref.  7  ]  for  transition  over  separation  bubbles.  They  used  for  Gv  the  ex- 
pression: 

c   =         lQQ 

y         exp(0.997V)    ' 

Where  Tu  is  the  turbulence  level  in  the  free  flow.  This  intermittency  model  was  also  in- 
vestigated durins  the  work  on  this  thesis. 


5.     Transition 

The  prediction  of  transition  from  laminar  to  turbulent  flow  is  very  difficult  and 
has  to  rely  on  empirical  correlations.  The  relation  used  here  to  predict  the  onset  of 
transition  is  a  combination  of  Michel's  method  and  the  e9  method,  and  is  given  by 
Cebeci  and  Bradshaw  [Ref.  2  ,  p.    153]: 

r^    =  1.174(1  +   4^  V46 


1  tt 


cirr 


Where: 

1.  R6ir  is  the  Reynolds  number  based  on  the  momentum  thickness  at  the  onset  of 
transition. 

2.  R,    is  the  Revnolds  number  based  on  x  at  the  onset  of  transition. 

In  the  computer  code,  if  a  laminar  separation  is  detected  before  transition  oc- 
curs, the  onset  of  transition  is  assumed  at  the  point  of  laminar  separation. 
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III.     DESCRIPTION  OF  THE  COMPUTER  CODE 

The  computer  code  used  here  to  investigate  cascade  flows  was  written  by  Cebeci, 
and  is  based  on  the  numerical  formulation  that  was  outlined  in  the  previous  chapters. 
In  this  chapter  the  general  structure  and  the  major  subroutines  of  the  code  will  be  de- 
scribed. 

A.  GENERAL  STRUCTURE  OF  THE  MAIN  PROGRAM 

The  main  program  reads  in  the  cascade  data  (blade  coordinates,  spacing  and  stagger 
angle),  the  flow  data  (inlet  angle  and  Reynolds  number),  and  transition  parameters.  The 
transition  onset  on  each  surface  of  the  blade  can  be  computed  by  the  program,  or  can 
be  input  by  the  user.  The  intermittency  parameter  G  should  be  specified  by  the  user. 

The  program  then  calls  subroutine  POTNL  to  compute  the  outer  inviscid  flow  field 
for  the  first  cycle.  The  output  of  subroutine  POTNL  is  the  external  velocity  distribution 
on  the  surface  of  the  blades.  This  velocity  distribution  is  then  transferred  to  subroutine 
CASBLP,  which  calculate  the  boundary  layer  flow. 

Subroutine  CASBLP  returns  the  displacement  thickness  distribution  and  the  blow- 
ing velocity  distribution  on  the  blades  to  the  main  program.  This  data  is  then  transferred 
back  to  subroutine  POTNL  to  the  next  cycle  of  calculations. 

The  program  repeats  the  cycles  of  calculations  by  calling  the  two  subroutines,  until 
the  specified  number  of  cycles  is  reached,  or  until  a  convergence  criterion  is  satisfied. 

B.  DESCRIPTION  OF  THE  SUBROUTINES 
I.     Subroutine  POTNL 

This  subroutine  solves  the  inviscid  outer  flow  by  using  the  panel  method.  The 
subroutine  calculates  the  influence  coefficients  and  calculates  the  velocities  subject  to  the 
boundary  conditions. 

The  velocities  are  evaluated  on  the  displaced  surface  (the  surface  created  by 
adding  the  displacement  thickness  to  the  original  surface  of  the  blade).  The  input  to  this 
subroutine  includes  the  cascade  geometry,  the  blowing  velocity  and  the  displacement 
thickness  (for  the  first  cycle  both  the  displacement  thickness  and  the  blowing  velocity 
are  taken  to  be  zero). 
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2.  Subroutine  CASBLP 

This  subroutine,  called  by  the  MAIN  program,  receives  the  blade  geometry  and 
the  velocity  distribution  as  input. 

It  transforms  the  x.y  blade  coordinates  to  the  chordwise  tangential  coordinates 
and  smooths  the  velocity  data  (during  the  work,  on  this  thesis  it  was  found  that 
smoothing  the  velocity  data  prevents  the  detection  of  the  separation  bubble  near  the 
leading  edge,  and  therefore  it  was  eliminated).  The  subroutine  then  calls  subroutine 
COMPBL  for  further  calculations. 

3.  Subroutine  COMPBL 

This  subroutine  finds  the  stagnation  point  and  controls  the  generation  of  the 
boundary  layer  calculation  grid  for  each  surface  (the  grid  starts  at  the  stagnation  point 
and  includes  91  points  in  the  chordwise  direction  for  the  upper  surface  and  71  points  on 
the  lower  surface). 

The  subroutine  then  calls  subroutine  BL2D  which  calculates  the  boundary  layer 
parameters  for  each  surface  (BL2D  is  called  twice,  first  for  the  upper  surface  and  then 
for  the  lower  surface). 

4.  Subroutine  BL2D 

This  subroutine  computes  the  displacement  thickness  and  the  blowing  velocity 
and  returns  them  back  to  the  calling  subroutine  (COMBL)  in  arrays  compatible  with  the 
potential  flow  calculations  (one  array  that  contains  all  the  points  of  the  blade,  first  the 
lower  surface  starting  at  the  trailing  edge  and  proceeding  forward,  and  then  the  upper 
surface,  starting  at  the  leading  edge  and  proceeding  backwards). 

BL2D  calls  the  following  subroutines: 

1.  Subroutine  INPUT  which  calculates  the  following: 

a.  NS,  the  switching  point  between  direct  and  interactive  boundary  layer  calcu- 
lations (this  point  is  set  at  the  first  pressure  peak  when  the  blade  is  scanned  from 
leading  edge  towards  the  trailing  edge) 

b.  NTR,  transition  location  (only  if  the  transition  location  is  an  input.  Otherwise 
it  is  calculated  by  subroutine  TRNS). 

c.  GMTR,  the  distribution  of  the  intermittency  factor  ytr. 

In  addition  this  subroutine  generates  the  boundary  layer  grid  in  the  r\  direction  and 
the  initial  velocity  profile,  by  calling  subroutine  INTL. 

2.  CALCIJ,  calculates  the  c„  coefficients  used  in  the  Hilbert  integral  approximation. 

3.  EDDY,  calculates  the  eddy  viscosity  (called  only  after  transition  has  been  de- 
tected). 
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4.  COEFTR.    calculates    the   coefficients    of  the    boundary    layer   finite    difference 
equations  in  transformed  form  (  for  the  direct  method  calculations). 

5.  SOLVE3,  solves  the  linearized  boundary  layer  equations  for  the  F,U  and  V  vari- 
ables by  computing  the  increments  <5F,  <5U  and  <5V 

The  subroutine  then  checks  the  convergence  of  the  Newton  iterations  and  re- 
peats the  calculations  if  needed.  If  the  subroutine  detects  flow  separation  or  if  it  reaches 
the  switching  point  NS,  subroutine  MAIN2  is  called  for  the  interactive  method  calcu- 
lations. Otherwise,  the  subroutine  proceeds  to  the  next  chordwise  point  of  the  grid  (XX) 
and  repeats  the  calculations. 
5.     Subroutine  MAIN2 

This  subroutine  calculates  the  boundary  layer  parameters  by  the  interactive 
method.  The  subroutine  performs  the  following  steps: 

1.  It  first  calls  subroutines  JOINT  and  COMGI  to  compute  the  interaction  coeffi- 
cients. 

2.  In  regions  of  laminar  flow  it  calls  the  following  subroutines: 

a.  COEF.  which  calculates  the  coefficients  of  the  boundary  layer  finite  differences 
equations. 

b.  SOLV4,  solves  for  the  variables  F,  U,  V  and  W  by  computing  the  increments 
<5F,  <5U  ,  <5V  and  <5W  . 

c.  TRANS,  to  check  if  the  condition  for  transition  is  satisfied  (it  also  checks  for 
laminar  separation  and  initiates  transition  at  the  point  of  laminar  separation  if 
it  is  detected). 

The  subroutine  then  checks  for  convergence  of  Newton  iterations  and  repeats  the 
calculations  as  needed. 

3.  In  regions  of  turbulent  flows  the  subroutine  calls  the  following  subroutines: 

a.  EDDY,  to  compute  the  eddy  viscosity  parameter  B   (B  =  1  +  vr/v  ). 

b.  COEF  and  SOLV4,  the  same  as  for  laminar  flow. 


6.     Subroutine  OUTPUT 

This  subroutine  computes  the  boundary  layer  parameters.  It  is  called  with  a 
parameter  "INDEX"  which  determines  the  type  of  calculations: 

1.   For   INDEX  =1    the   computations   relates   to   transformed    coordinates   (direct 
boundary  layer  method)  using  the  relations: 

2V(1.2)B(1,2) 
c/=    = • 


■JKX 
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Vw  -  uJ^-V(\,2), 


D  =  ueS  jRe  , 


.  ♦ 


-g=r{t,{NF)-F(NP))t 


NP 


6    = 


7=2 

Where  V(l,2)  and  B(  1,2)  are  the  velocity  gradient  and  the  viscosity  parameter  at  the 
surface,  respectively,  and  ?/(NP)  and  F(NP)  are  r\  and  F  evaluated  at  the  edge  of  the 
boundary  layer. 

For  INDEX  =  2  the  subroutine  calculates  the  boundary  layer  parameters  for  semi- 
-transformed  coordinates  (interactive  boundary  layer  method)  using  the  relations: 


c/ 

2V(1,2)B(1,2) 

JJFjK,  [W(A7>)]2 

V(1.2) 

si* 

ue  =  U(AT) , 

D  =  (Uj  -  F)v7  , 

<5*  =  hi-—JJx  . 

For  NX  >  NTR  (after  transition  has  been  detected),  subroutine  SMPSON  is 
called  (subroutine  SMPSON  calculates  the  coefficient  of  the  outer  region  eddy  viscosity). 
The  subroutine  then  prints  out  the  velocity  profiles  at  the  required  stations. 

7.  Subroutine  TRANS 

This  subroutine  calculates  the  transition  location  based  on  the  Michel  criterion 
or  based  on  laminar  separation  (whichever  occurs  first).  If  transition  has  been  detected 
the  intermittency  distribution  is  calculated  for  all  the  remaining  points  of  the  surface. 

8.  Subroutine  FILLUP 

This  subroutine  increases  the  number  of  points  in  the  boundary  layer  grid  (in 
the  r]  direction)  as  needed.  It  also  fills  up  the  arrays  of  F,  U,  B,  W  and  V  between  the 
edge  of  the  boundary  layer  to  the  end  of  the  arrays  (with  V=0,  \V,B  and  U,  with  the  last 
values  that  they  had  in  the  edge  of  the  boundary  layer  and  F  as  the  integral  of  U). 


16 


9.  Subroutine  EDDY 

This  subroutine  calculates  the  eddy  viscosity  using  the  Cebeci-Smith  two  layer 
eddy  viscosity  formula.  It  receives  the  vectors  U.V  and  r\  at  a  point  and  computes  the 
viscosity  vector  B. 

10.  Subroutine  INTL 

This  subroutine  generates  the  boundary  layer  grid  in  the  r\  direction.  It  sets  the 
number  of  grid  points  and  generates  the  initial  velocity  profile. 
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IV.     RESULTS  AND  DISCUSSION 

The  viscous  inviscid  interaction  code  was  run  with  several  cascades  on  which  ex- 
perimental data  is  available.  In  order  to  enable  a  thorough  comparison  between  exper- 
imental results  and  the  computed  results,  a  very  detailed  experimental  data  base  is 
needed.  The  data  should  include  measurements  of  the  boundary  layer  development  along 
the  blade,  velocity  profiles  along  the  boundary  layer,  transition  location  and  distribution, 
flow  separation,  and  external  velocity  distribution. 

Unfortunately,  very  few  cascade  experiments  has  been  performed,  which  obtained 
the  required  data  with  sufficient  accuracy,  due  mostly  to  the  lack  of  appropriate  meas- 
urement equipment.  Only  recently,  with  the  introduction  of  non-interfering  methods 
like  the  Laser  Doppler  Velocimeter  (LDV),  the  required  data  can  be  measured  accu- 
rately. 

Recently  an  experiment  involving  the  investigation  of  a  linear  compressor  cascade 
of  Controlled  Diffusion  Blading  (which  will  be  referred  here  as  the  CD  cascade)  has  been 
carried  out  by  Elazar  [Ref.  8].  Most  of  the  work  in  the  present  thesis,  involves  compar- 
ison of  the  computer  code  results  with  Elazar's  experimental  results. 

Other  cascades  that  were  investiaated  are: 

1.  A  shockless,  supercritical  airfoil  cascade,  designed  in  1974  by  Korn  in  cooperation 
with  Pratt  &  Whitney  Aircraft  (referred  here  as  the  P  &,  W  cascade).  The  exper- 
imental results  of  the  cascade  were  obtained  from  a  report  by  Hobbs,  Wagner. 
Dannenhoffer  and  Dring  [Ref.  9]. 

2.  Stator  blade  of  a  single  stage  axial  compressor  (referred  here  as  the  C4  cascade). 
The  blade  profile  is  the  British  C4  section  (10%  thickness)  on  a  circular  arc  camber 
line.  The  experiment  has  been  performed  by  Walker  [Ref.  10].  The  detailed 
boundary  layer  measurements  are  not  presented  in  the  report  and  were  obtained 
directly  from  the  author. 

The  code  failed  to  run  with  two  other  cascades: 

1.  A  highly  loaded,  double  circular  arc  blade  with  a  sharp  leading  edge  and  a  sharp 
trailing  edse,  used  in  a  compressor  cascade  that  was  investigated  bv  Deutsch  and 
Zierk[Ref.  11]. 

2.  V2  double  circular  arc  blade,  highly  loaded  cascade.  This  cascade  was  investigated 
by  Hoheisel  and  Seyb  [Ref.  12]. 

In  both  cases  the  code  calculated  the  potential  flow  successfully  but  failed  in  trying 
to  compute  the  first  cycle  of  the  boundary  layer  calculations. 
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A.     CD  CASCADE 

The  experimental  data  for  the  CD  cascade  was  obtained  at  M  =  0.25,  Rf  =700000 
and  at  three  inlet  angles:  40°  (the  design  condition).  43.4°  and  46°  .  The  spacing  was 
0.6  of  the  chord  and  the  stagger  angle  14.27°  .  A  general  layout  of  the  cascade  i^  shown 
in  Figure   1  on  page  20. 

The  following  observations  were  concluded  from  the  experiment: 

1.  A  separation  bubble  exists  near  the  leading  edge  on  the  upper  surface  at  all  the  inlet 
angles.  The  bubble  became  larger  at  increased  inlet  angles. 

2.  Transition  from  laminar  to  turbulent  flow  occurred  above  the  separation  bubble 
(on  the  upper  surface). 

3.  Transition  on  the  lower  surface  occurred  at  midchord. 

4.  The  boundary  layer  thickness  on  the  upper  surface  increased  with  inlet  angle,  and 
reached  a  thickness  of  15%  chord  at  the  highest  inlet  angle.  The  boundary  layer 
thickness  on  the  lower  surface  did  not  change  significantly  with  inlet  angle. 

5.  The  turbulent  boundary  layer  on  both  surfaces  remained  fully  attached  at  all  the 

inlet  angles. 

1.     Transition  location  and  intermittency  distribution 

The  effects  of  the  transition  location  and  the  intermittency  factor  were  investi- 
gated. The  code  was  first  run  with  the  transition  location  calculated  bv  the  code,  and 
with  several  values  of  the  intermittency  factor  Gy.  It  was  found  that  the  code  did  not 
run  with  Gy  =  1200  (which  is  the  value  used  usually  for  high  Reynolds  numbers).  The 
highest  value  of  Gv  with  which  the  code  run  successfully  was  900. 

The  code  failed  to  predict  the  separation  bubble  on  the  upper  surface,  and  pre- 
dicted laminar  separation  at  7S%  chord  on  the  lower  surface  (which  did  not  occur  in  the 
experiment).  Transition  on  the  upper  surface  occurred  at  41%  chord  (detected  by 
Michel's  criterion)  and  at  78%  chord  on  the  lower  surface  (at  laminar  separation). 

The  shape  factor  computed  by  the  code  was  compared  to  the  experimental  re- 
sults. As  can  be  seen  in  Figure  2  on  page  21  the  shape  factor  as  predicted  by  the  code 
deviates  substantially  from  the  actual  results,  due  mainly  to  the  different  transition  lo- 
cation. 

On  the  lower  surface,  as  can  be  seen  in  Figure  3  on  page  22  the  shape  factor 
deviates  even  more  from  the  experimental  results.  In  this  figure  the  effect  of  changing 
the  intermittency  factor  G,  can  be  seen.  For  both  the  extreme  values  of  Gv  .  10  and  900, 
the  computed  shape  factor  curve  is  far  from  agreement  with  the  actual  results. 
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Figure   1.      Controlled  Diffusion  cascade 
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Figure  2.      Shape  factor  comparison  on  the  tipper  surface:     Transition  computed  by 
the  code  (0  =  40°) 
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Figure  3.      Shape  factor  comparison  on  the  loner  surface:     Transition  computed  by 
the  code  (/?  =  40°). 
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Since  the  code  did  not  predict  the  existence  of  the  separation  bubble  on  the 
upper  surface,  it  was  decided  to  try  to  eliminate  the  smoothing  of  the  external  velocity 
as  computed  by  the  potential  flow  subroutine  (originally,  the  velocities  were  smoothed 
in  subroutine  CASBLP  prior  to  boundary  layer  calculations).  It  was  found  that  without 
smoothing  the  velocities  a  small  separation  bubble  is  predicted  by  the  code  at  4%  of 
chord.  The  onset  of  transition  is  set  by  the  code  at  the  beginning  of  the  separation 
bubble. 

The  code  was  run  with  unsmoothed  velocities  with  two  values  of  Gy,  10  and  900. 
The  shape  factor  behavior  can  be  seen  in  Figure  4  on  page  24.  Changing  the  value  of 
Gy  did  not  change  the  shape  of  the  curve  much,  and  generally  the  shapes  of  the  com- 
puted and  the  experimental  curves  look  alike. 

The  elimination  of  the  velocity  smoothing  in  the  code,  also  affects  the  thickness 
of  the  boundary  layer.  In  Figure  5  on  page  25  the  displacement  thickness  is  plotted  for 
both  cases  (with  and  without  the  velocity  smoothed).  Without  smoothing,  the  displace- 
ment thickness  is  much  thicker,  especially  on  the  rear  half  of  the  blade,  which  is  closer 
to  the  actual  results. 

The  effect  of  changing  the  intermittency  distribution  to  the  one  used  by  Rodi 
and  Schonung  [Ref.  7]  was  investigated.  It  was  found,  as  can  be  seen  in  Figure  6  on  page 
26  that  the  effect  of  the  new  model  is  equivalent  to  using  Gv  in  the  present  model. 

On  the  lower  surface  it  was  necessary  to  run  the  code  with  transition  as  input, 
to  get  reasonable  results,  as  can  be  seen  in  Figure  7  on  page  27  for  transition  input  at 
21%  of  chord. 

At  the  of!  design  conditions  (inlet  angles  of  43.4°  and  46°)  a  similar  behavior  of 
the  transition  has  been  observed,  as  can  be  seen  for  example  in  Figure  8  on  page  28  for 
the  upper  surface  and  in  Figure  9  on  page  29  for  the  lower  surface,  both  at  /)  =  46°. 
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Figure  4.      Shape  factor  on  the  upper  surface  without  velocity  smoothing. 
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Figure  5.      Displacement  thickness:     The  effect  of  velocity  smoothing. 
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Figure  6.      The  effect  of  the  intermittency  model:     Upper  surface,  f]  =  40c 
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Figure  7.      Shape  factor  on  the  lower  surface  with  transition  input  at  21%. 
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Figure  8.      Shape  factor  at  /J  =  -46°  on  the  upper  surface 
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Figure  9.      Shape  factor  at  /I  =  46°  on  the  loner  surface. 
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2.     External  Velocity  Distribution 

The  external  velocity  distribution,  computed  by  the  code  using  the  interaction 
law,  was  compared  to  the  experimentally  measured  velocities.  It  was  found  that  in  gen- 
eral, the  velocities  measured  experimentally,  were  higher  than  those  computed  by  the 
code  for  all  the  inlet  angles. 

There  are  two  possible  sources  to  the  discrepancy  in  the  velocities: 

1.  The  computer  code  calculates  pure  2--D  flows.  In  the  experiment  the  flow  was 
observed  to  accelerate  due  to  the  effect  of  the  boundary  layer  on  the  side  walls  (a 
3-D  effect).  This  effect  was  calculated  in  the  experiment  and  is  referred  to  as  the 
AVDR  correction  [Ref.  8,  p.43]. 

2.  The  flow  accelerates  due  to  the  thickening  of  the  boundary  layer.  Since  the 
boundary  layer  as  computed  by  the  code  is  substantially  thinner  than  the  actual 
boundary  layer  (as  will  be  discussed  in  the  next  section)  the  external  velocities 
predicted  by  the  code  are  smaller. 

To  compensate  for  the  first  error  source,  all  the  computed  velocities  were  com- 
pared to  the  experimental  velocities  corrected  by  the  AVDR  correction.  The  comparison 
between  the  velocities  can  be  seen  in  Figure  10  on  page  31  for  ft  =  40°,  in  Figure  11 
on  page  32  for  /?  =  43.4°  and  in  Figure  12  on  page  33  for  /?  =  46°. 

It  can  be  seen  from  the  figures  that  the  difference  between  the  computed  and  the 
experimental  velocities  is  larger  on  the  lower  surface.  The  reason  might  be  the  method 
with  which  the  correction  to  the  inviscid  velocity  is  computed.  The  assumption  on  which 
the  interaction  law  is  based,  is  that  only  sources  (representing  the  viscous  effects)  on  the 
surface  being  considered,  affect  the  local  velocity.  In  reality,  the  boundary  layer  on  both 
surfaces  affects  the  local  velocity  (because  the  boundary  layer  developed  on  the  upper 
surface  of  a  blade,  causes  a  velocitv  disturbance  that  is  felt  on  the  lower  surface  of  the 
adjacent  blade). 

Since  the  boundary  layer  on  the  lower  surface  is  much  thinner,  its  effect  on  the 
velocity  on  the  upper  surface  is  much  smaller  than  the  effect  of  the  upper  surface 
boundarv  laver  on  the  lower  surface  velocitv. 
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Figure   10.      External  velocity  at  ft  =  40s 
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Figure   II.      External  velocity  at /I  =  43. 4C 
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Figure   12.      External  velocity  at  /?  =  46c 
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3.     Boundary  Layer  Thickness 

The  boundary  layer  thickness  as  computed  by  the  code  was  compared  with  the 
experimental  results  by  comparing  the  displacement  thicknesses. 

It  was  found  that  on  the  lower  surface  the  computed  and  the  actual  displace- 
ment thickness  agree  quite  well,  as  can  be  seen  in  Figure  13  on  page  35  for  /?  =  40°, 
Figure  14  on  page  36  for  /?  =  43.4°  and  in  Figure  15  on  page  37  for  /?  =  46°  . 

On  the  upper  surface  the  displacement  thicknesses  computed  by  the  program 
are  significantly  thinner  than  those  measured  experimentally.  The  difference  between  the 
computed  and  the  actual  thickness  increases  along  the  blade  and  it  increases  with  in- 
creased inlet  angle.  It  was  found  that  by  using  a  different  expression  for  the  inner  region 
eddy  viscosity  (as  mentioned  in  chapter  II),  the  displacement  thickness  can  be  increased, 
but  the  difference  between  the  actual  and  the  computed  thickness  is  still  substantial,  es- 
pecially at  the  higher  inlet  angles.  Figure  16  on  page  38,  Figure  17  on  page  39  and 
Figure  18  on  page  40  shows  the  displacement  thickness  on  the  upper  surface  for  the 
three  inlet  aneles,  with  the  original  and  the  modified  eddv  viscosity  models. 

The  large  error  in  the  prediction  of  the  boundary  layer  thickness,  can  be  the  re- 
sult of  several  reasons: 

1.  The  transition  model  used  in  the  code,  sets  the  onset  of  transition  at  the  first  point 
of  laminar  separation.  It  causes  rapid  transition  to  turbulent  flow  which  reattaches 
immediately,  resulting  in  a  very  small  separation  bubble  compared  to  the  bubble 
observed  in  the  experiment. 

2.  The  turbulent  model  used  in  the  code  could  be  inaccurate.  It  was  derived  based  on 
empirical  data  obtained  in  single  airfoil  experiments  and  not  with  cascades.  In  ad- 
dition the  present  model  does  not  include  the  effects  of  the  free  stream  turbulence 
(that  was  relatively  high  in  the  experiment). 

3.  The  boundary  layer  as  measured  in  the  experiment  was  quite  thick,  especially  at  the 
higher  inlet  angles  (it  reached  15%  of  the  chord  at  /?  =  46°).  Such  a  thick  bound- 
ary layer  may  violate  the  basic  assumptions  on  which  the  boundary  layer 
equations,  and  the  interaction  law,  were  based  (especially  when  the  spacing  be- 
tween the  blades  is  small,   60%  chord  in  this  case). 

It  was  suggested  that  one  of  the  possible  reasons  to  the  inaccurate  prediction 
of  the  boundary  layer  is  the  blunt  trailing  edge  of  the  blade,  that  might  cause  difficulties 
in  the  computations.  A  modified  blade,  with  a  sharp  trailing  edge  has  been  run,  and  the 
displacement  thickness  distribution  can  be  seen  in  Figure  19  on  page  41.  As  can  be  seen 
in  the  figure  the  sharp  trailing  edge  affects  only  the  boundary  layer  adjacent  to  the 
trailing  edge,  and  therefore  cannot  provide  an  explanation  to  the  difference  between  the 
actual  and  the  computed  displacement  thickness. 
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Figure   13.      Displacement  thickness  on  the  ltmer  surface  (/?  =  40°  ) 
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Figure   14. 


Displacement  thickness  on  the  loner  surface  (/J  =  43.4°  ) 
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Figure    15.      Displacement  thickness  on  the  hmer  surface  (//  =  46°  ) 
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Figure   16.      Displacement  thickness  on  the  upper  surface  (/J  =  40°  ) 
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Figure   17.      Displacement  thickness  on  the  upper  surface  (/?  =  43.4°  ) 
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Figure   18.      Displacement  thickness  on  the  upper  surface  (/J  =  46°  ) 
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Figure   19.      The  effect  of  sharp  trailing  edge. 
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4.     Comparison  to  a  Navier  Stokes  Code 

A  limited  comparison  of  the  experimental  and  the  computed  results  with  a 
Navier  Stokes  (N.S.)  code  calculations  has  been  performed.  The  N.S.  code  has  been  de- 
veloped and  run  by  S.  J.  Shamroth  of  Scientific  Research  Associates  inc.  in  cooperation 
with  Pratt  and  Whitney  Aircraft. 

Since  the  N.S.  code  does  not  compute  the  displacement  thickness,  the  velocity 
profiles  near  the  surface  of  the  blade  were  compared.  The  comparisons  were  made  at 
90%  chord  on  the  suction  surface  for  all  three  inlet  angles. 

At  the  design  point.  /?  =  40°,  shown  in  Figure  20  on  page  43.  both  the  inter- 
active code  and  the  N.S.  code  failed  to  predict  accurately  the  actual  velocity  profile.  In 
this  case  the  interactive  code  seems  to  yield  somewhat  better  results  than  the  N.S.  code. 

At  the  higher  inlet  angles.  /?  =  43.4°  and  /?  =  46°,  shown  in  Figure  21  on  page 
44  and  in  Figure  22  on  page  45  respectively,  the  N.S.  calculations  show  significantly 
better  agreement  with  the  experimental  results  than  the  interactive  code. 

From  these  comparisons,  it  can  be  seen  that  the  interactive  code  deviation  from 
the  actual  results  increases  with  increased  inlet  ancle  (increased  loadine  of  the  cascade), 
whereas  the  N.S.  code  deviation  seems  to  decrease  with  increased  inlet  ansle. 
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Figure  20.      The  results  of  the  N.  S.  code  at  fi  -  40° 
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Figure  21.      The  results  of  the  N.  S.  code  at  P  =  43. 4C 


44 


o 


CO 

d 


CO 

w  ° 


d 


0.00    0.02 


LEGEND 
□  EXPERIMENTAL 


I NTE RACTIVE  CODE 

N.S.CODE 


0.04  0.06  0.00  0.10 

D/C 

BETA=46 


0.12 


Figure  22.      The  results  of  the  N.  S.  code  at  p  =  46c 
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B.     P  &  W  CASCADE 

The  experimental  data  for  the  P  &  \V  cascade  was  obtained  at  inlet  flow  angle  of 
52°,  at  VI  =  0.11,  and  Reynolds  number  of  478000.  The  cascade  had  a  stagger  angle  of 
15.75°  and  0.7  spacing.  A  general  layout  of  the  cascade  is  shown  in  Figure  23  on  page 
47. 

A  comparison  of  the  computed  and  the  measured  pressure  coefficients  on  the  blade 
is  shown  in  Figure  24  on  page  48.  There  is  a  good  agreement  between  the  computed 
and  the  measured  CP. 

The  displacement  thickness  was  measured  in  the  experiment  only  at  96.8%  of  chord. 
This  measurement  is  compared  to  the  computed  results  in  Figure  25  on  page  49.  As  can 
be  seen,  the  computed  and  the  measured  data  agree  almost  perfectly  on  the  lower  sur- 
face, and  quite  well  on  the  upper  surface.  The  difference  observed  on  the  upper  surface 
is  caused  by  the  early  prediction,  by  the  code,  of  trailing  edge  separation,  a  short  dis- 
tance upstream  of  the  actual  location.  This  can  also  be  observed  when  comparing  the 
velocity  profiles  at  that  point,  in  Figure  26  on  page  50.  The  computed  velocity  curve 
shows  a  small  zone  of  reversed  flow  near  the  surface  of  the  blade.  This  reversed  flow 
could  be  the  result  of  a  too  early  prediction  of  trailing  edge  separation  by  the  code,  or 
it  could  have  existed  in  the  actual  flow  but  not  detected  because  of  its  size. 
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Figure  23.      Pratt  &  Whiteny  cascade 
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Figure  24.      Comparison  of  pressure  coefficient. 
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Figure  25.       Displacement   thickness   comparison:     Experimental   data    shown   at 
96.8%  chord. 
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Figure  26.      Velocity  profile  at  96.8%  chord  on  the  upper  surface. 
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C.     C4  CASCADE 

The  C4  cascade  has  a  stagger  angle  of  29.5°,  a  camber  angle  of  31.1°  and  spacing 
of  0.992.  It  has  been  tested  at  Reynolds  numbers  of  about  200000,  and  inlet  angles  of 
34.1°  to  47.7°  {which  corresponds  to  incidence  angles  of  -10.9°  to  2.7C).  The  general 
layout  of  the  cascade  is  shown  in  Figure  27  on  page  52.  A  computer  code  that  generates 
the  coordinates  of  the  blade  and  a  summary  of  some  experimental  results  are  given  in 
Appendix  B. 

The  code  was  run  with  the  intermittency  constant  Gv  =  10.  Higher  values  of  Gv 
(above  100)  caused  numerical  problems  in  the  code.  The  onset  of  transition  was  first 
taken  at  the  point  where  it  was  observed  in  the  experiment.  At  the  lower  inlet  angle  it 
seems  that  a  better  agreement  with  the  experimental  results  can  be  obtained  by  delaying 
the  onset  of  transition  but  trying  to  implement  it  resulted  in  numerical  breakdown  of  the 
computation.  At  the  higher  inlet  angles,  better  agreement  with  the  experimental  results 
was  achieved  by  initiating  the  transition  earlier  (at  26%  chord  for  p  =  45. 6'  and  at  21% 
for  /?  =47.7°  as  compared  to  44%  and  36%  chord  as  observed  in  the  experiment). 
1.     Displacement  Thickness 

Comparisons  of  the  experimental  data  to  the  computed  displacement  thickness 
are  shown  in  Figure  28  on  page  53  for  inlet  angle  of  34.1°,  in  Figure  29  on  page  54  for 
inlet  angle  of  36.3°,  in  Figure  30  on  page  55  for  inlet  angle  of  45.6°  and  in  Figure  31 
on  page  56  for  inlet  angle  of  47.7°. 

As  can  be  seen  in  the  figures,  there  is  a  good  agreement  between  the  actual  and 
the  computed  results  at  the  two  lower  angles  (ft  =  34.1°  and  /?  =  36.3°,  in  which  the 
incidence  angles  were  negative).  At  the  two  higher  angles.  /?  =  45.6°  and/)  =  47.7°  the 
computed  results  agree  with  the  actual  results  up  to  about  70%  chord,  and  then  the 
displacement  thickness  predicted  by  the  code  becomes  much  thicker  than  the  actual  one. 

The  code  predicted  a  large  flow  separation  area  starting  at  about  70%  chord  at 
the  lower  inlet  angles,  and  at  about  46%  chord  at  the  higher  inlet  angles.  This  flow 
separation  was  not  obsc-ved  in  the  experiment.  The  discrepancies  between  the  computed 
and  the  actual  results  behind  60%  to  70%  chord  can  be  explained  by  the  inaccurate 
calculations  by  the  code  due  to  the  large  separated  areas.  When  the  code  encounters 
separation,  several  approximations  are  made  (like  the  FLARE  approximation)  based  on 
the  assumption  that  the  separated  area  is  small.  When  the  separated  area  is  large,  these 
approximations  may  result  in  inaccurate  prediction  of  the  flow  field. 
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Figure  27.      C4  Cascade 
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Figure  28.       C4  cascade  at  /J  =   3-4.1°:     Displacement  thickness  comparison  with 
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Figure  30.       C4  cascade  at  /?  =  45.6°:     Displacement  thickness  comparison  with 
computed  results  (G=  10). 
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Figure  31.      C4  cascade  at  /?  =   47.7°:     Displacement  thickness  comparison  with 
computed  results  (G=  10). 
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2.     External  Velocity  and  Velocity  Profiles  Comparisons 

A  comparison  of  the  external  velocity  on  the  upper  surface  of  the  blade  is  shown 
in  Figure  32  on  page  58  for  inlet  angle  of  45.6°  and  in  Figure  33  on  page  59  for  inlet 
angle  of  47.7°.  It  can  be  seen  that  there  is  a  good  agreement  between  the  experimental 
and  the  computed  results  up  to  about  80%  chord.  Near  the  trailing  edge  the  computed 
results  deviate  from  the  experimental  results  due  to  the  inaccuracy  in  the  calculations 
of  the  displacement  thickness. 

A  comparison  of  the  velocity  profiles  in  the  boundary  layer  at  50%  chord  is 
shown  in  Fisure  34  on  page  60  for  inlet  angle  of  34.1°  and  in  Figure  35  on  page  61  for 
inlet  angle  of  36.3°.  The  agreement  between  the  calculated  velocity  profiles  and  the 
measured  velocity  profiles  is  very  good. 
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Figure  32.      C4  cascade  at  (l  =  -15.6°:     External  velocity  distribution. 
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Figure  33.      C4  cascade  at  fl  =  ALT:     External  velocity  distribution. 
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Figure  34.      C4  cascnde  at  /?  =  34.1°:     Velocity  profile  at  50%  chord. 
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V.     CONCLUSIONS  AND  RECOMMENDATIONS 

A.     CONCLUSIONS 

The  interactive  viscous  inviscid  computer  code,  has  been  investigated  by  comparing 
its  predictions  of  boundan.  layer  parameters  to  experimental  data. 

It  has  been  found  that  the  code  yields  reasonable  results  for  lightly  loaded  cascades. 
but  the  prediction  of  the  boundary  layer  thickness  on  the  suction  surface  of  highly 
loaded  cascades  deviates  significantly  from  experimentally  measured  data.  In  two  cases 
involving  highly  cambered  cascade  blades,  with  sharp  leading  edge,  the  code  failed  to 
run.  It  was  also  found  that  the  prediction  of  external  velocity  distribution  on  highly 
loaded  cascades  was  inaccurate. 

The  main  reasons  to  the  discrepancies  in  the  prediction  of  the  boundary  layer 
thickness  seem  to  be: 

1.  Inaccuracy  in  predicting  flow  parameters  in  regions  of  large  flow  separation  (due 
to  inadequate  transition  model  and  approximations  made  in  calculating  the  How  in 
separated  areas). 

2.  Inaccurate  turbulence  modelling. 

3.  Possible  violation  of  the  basic  assumptions  of  the  boundary  layer  theory  in  areas 

oi  very  thick  boundary  layers. 

4.  The  wake  is  not  calculated  by  the  code.  The  result  is  inaccurate  flow  prediction 
near  the  trailing  edge. 

The  inaccuracy  in  the  prediction  of  the  external  velocity  distribution  in  highly  loaded 
cascades  is  due  to  the  interaction  law.  which  does  not  account  for  the  presence  of  adja- 
cent blades. 


B.     RECOMMENDATIONS 

The  recommended  steps  in  order  to  improve  the  code  are: 

1.  Improving  the  interaction  law  by  assuming  a  distribution  of  sources  on  the  actual 
surface  (instead  of  the  assumption  of  a  Hat  plate),  letting  the  correction  term  to  the 
external  velocity  vary  across  the  boundary  layer  and  distributing  sources  on  the 
adjacent  blade  as  well  for  better  modelling  of  the  boundary  layer  effect  on  the  ex- 
ternal velocitv. 


2.  Changes  to  the  derivation  of  the  boundary  layer  equations  should  be  investigated 
to  allow  a  better  treatment  of  thick  boundary  layers  (like  omitting  the  assumption 
ofcpjcy  =  0  across  the  boundary  layer). 

3.  Different  turbulence  models  should  be  investigated. 

4.  The  wake  should  be  included  in  the  calculations. 
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APPENDIX  A.     COMPUTER  CODE  LISTING 


C** 

c** 
c** 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


VISCOUS-INVISCID  INTERACTION  PROGRAM  FOR  CASCADE  FLOWS 


•vWcsWo'c 


VERSION 
JANUARY  87 


3.  A 


THIS  VISCOUS-INVISCID  INTERACTION  METHOD,  CAPABLE  OF  COMPUTING  BOTH 
SINGLE  AIRFOIL  AND  CASCADE  FLOWS,  WAS  DEVELOPED  BY  CEBECI  AND 
COLLABORATEURS  AT  LONG  BEACH  STATE  AND  DOUGLAS  AIRCRAFT  COMPANY. 
THE  CODE  APPLIES  TO  INCOMPRESSIBLE,  2-DIMENSIONAL,  STEADY  FLOWS 
PAST  LINEAR,  ARBITRARILY  STAGGERED  CASCADES.  THE  METHODS  BASIC 
INGREDIENTS  INC'.UDE 

1.  A  FIRST  ORDL1  PANEL  METHOD  TO  SOLVE  LAPLACE'S  EQUATION, 

2.  A  FINITE  DIFFERENCE  SCHEME  TO  SOLVE  THE  BOUNDARY  LAYER  EQUATIONS 
SUBJECT  TO  DIRECT  OR  INTERACTIVE  BOUNDARY  CONDITIONS, 


A  STRONG 


FRACTION  MODEL  TO  COUPLE  VISCOUS  AND  INVISCID  FLOW 


RESULTS,  AND 
4.   A  ZERO  EQUATION,  ALGEBRAIC  TURBULENCE  MODEL  TO  ESTIMATE 
TURBULENT  SHEAR  STRESSES. 

IN  SUMMARY,  THE  CODE  WILL  PROVIDE,  FOR  ATTACHED  AS  WELL  AS  MODERATE- 
LY SEPARATED  FLOWS  PAST  SINGLE  AIRFOILS  OR  CASCADES,  THE  FOLLOWING 

1.  INVISCID  AND  VISCOUS  PRESSURE  DISTRIBUTIONS, 

2.  DISTRIBUTIONS  OF 

A.  LOCAL  SKIN  FRICTION  COEFFICIENT, 

B.  DISPLACEMENT  AND  MOMENTUM  THICKNESS,  AND 

3.  VELOCITY  PROFILES  ACROSS  THE  BOUNDARY  LAYER. 

MODIFICATIONS  SINCE  VERSION  3.0: 

1.  PRECISE  ASSIGNMENT  OF  BEGIN  OF  TRANSITION. 

2.  CORRECTION  OF  AN  ERROR  IN  THE  CALCULATION  OF  MOMENTUM  THICKNESS. 

3.  ADDITIONAL  PRINT  OPTION:  IP=-2  WILL  PROVIDE  AN  INPUT  FILE  (UNIT 
NUMBER  12)  FOR  THE  PLOTTING  ROUTINE. 


COMMON  /BLCO/  NX,NXT,NP,NPT,NTR, IT, INVRS ,NS, IP 

COMMON/ B LOW/ VN( 100) 

COMMON/BLIN/  TITLE(20) ,XC( 100) ,YC( 100) , ISG( 100) ,DELS( 100) , 
+  XCTR,XTR,ISTRP,ICYCLE,ICYTL,XCTRS(2),TRFIND(2) 

COMMON/CASCDE/ INLET , SP , S INGLE , ALPHAA , ALPHAI , STAG 

COMMON/TRN/  PGAMTR , OMEGA , RTHETB , RTRANB 

COMMON/PLOT/NVP( 2 ) , NXVP( 20 , 2 ) , ICC 

DIMENSION  XO( 100) ,YO( 100) ,X( 100) ,Y( 100) ,VCOM( 100) ,DLS( 100) , 
+  XS( 100) ,YS( 100) ,XSTGR( 100) , YSTGR( 100) ,DBPP( 100) 

DIMENSION  CASEID(20),XCTRI(2),ITRI(2),NBL(2) 

LOGICAL  SINGLE,TRFIND 

TRFIND(1)=  .FALSE. 

TRFIND(2)=  .FALSE. 


INT0001 

INT0002 

INT0003 

INT0004 

INT0005 

INT0006 

INT0007 

INTO 003 

INT0009 

INT0010 

INT0011 

INT0012 

INT0013 

INT0014 

INTOOlff 

INT001I 

intooi; 

INT001I 

INTOOI! 

INT002( 

INT0021J 

INT00221 

INT0023} 

INT0024J 

INT0023 

INT0023 

INT0021 

INT002d 

INT0029 

INT003d 

INT003U 

INT0032] 

INT0033 

INT0034 

INT0033 

INT0036] 

INT003A 

INT003S 

INT0039I 

INT0040! 

INT004l| 

INT0042I 

INT0043 

INT0044 

INT0045I 

INT0043 

INT004A 

INT0048 

INT0049 
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1 

5 


C 

10 

20 
C 

25 

27 


30 

C 

C40 

40 

C 

45 

C 

50 


(NXVP(I,1),I=1,NVP(1)) 
(NXVP(I,2),I=1,NVP(2)) 


80 

82 


I CASE   =  0 

READ(5,5,END=999)  TITLE 

FORMAT(20A4) 

ICASE  =  ICASE  +  1 

REWIND  3 

READ  (5,10) 

FORM AT (IX) 

READ  (5,20)  ITRI(1),ITRI(2),IRST,ICYTL,IP 

F0RMAT(16I5) 

READ  (5,10) 
READ  (5,25)  INLET, I STAG, ALPHA I , STAG, SP,PGAMTR, OMEGA 
FORMAT(2I5,5F10. 0) 

READ  (5,27)RN,XCTRI(1),XCTRI(2),ALPHAA 
F0RMAT(4E10. 0) 
IF  (IP.EQ.  -2)  THEN 

READ  (5,20)  NVP(1),NVP(2) 
IF  (NVP(l).NE.  0)  READ  (5,20) 
IF  (NVP(2).NE. 0)  READ  (5,20) 
END  IF 

IF  (ICASE  . EQ.  1)  READ  (5,20)  N,NI 
I READ  =  1 
I BLOW  =  1 
SINGLE  =  .FALSE. 

IF  (SP  .  LE.  0.0)  SINGLE  =  .TRUE. 
N  =  N  -  1 
Nl=  N  +  1 

IF  (ICASE  .GT.  1)  THEN 
Nl  =  N1SAVE 
N  =  N  -  1 
GOTO  53 
END  IF 
IF  (IREAD  .EQ.  1)  GO  TO  40 

READ  (5,10) 
READ  (5,3C)  (XO(I),  YO(I),  1=1, N+l) 
FORMAT(2FI0. 0) 
GO  TO  50 

READ(5,10) 
READ(5,45)  (XO(I)  ,  1=1, N+l) 

READ(5,10) 
READ(5,45)  (YO(I)  ,  1=1, N+l) 
FORMAT(6F10.  0) 

CONTINUE 

IF  (IP.EQ.  -2)  THEN 

WRITE( 12,20)  N+1,NVP(1),NVP(2),90,70,INLET 

IF  (NVP(l).NE.O)  WRITE(12,20)  (NXVP( I , 1) , I=1,NVP( 1) ) 

IF  (NVP(2).NE.O)  WRITE(12,20)  (NXVP( I ,2) , 1=1 ,NVP( 2) ) 

IF  ( INLET. NE. 1)  WRITE(12,80)  RN,ALPHAA 

IF  (INLET.EQ. 1)  WRITE(12,80)  RN,ALPHAI 

WRITE(12,82)  (XO(I),I=l,N+l) 

W"RITE(12,82)  (YO(I),I=l,N+l) 

FORMAT(2E15.5) 

F0RMAT(8F10.  6) 
END  IF 


INT00500 
INT00510 
INT00520 
INT00530 
INT00540 
INT00550 
INT00560 
INT00570 
INT00580 
INT00590 
INT00600 
INT00610 
INT00620 
INT00630 
INT00640 
INT00650 
INT00660 
INT00670 
INT00680 
INT00690 
INT00700 
INT00710 
INT00720 
INT00730 
INT00740 
INT00750 
INT00760 
INT00770 
INT00780 
INT00790 
INT00800 
INT00810 
INT00820 
INT00830 
INT00840 
INT00850 
INT00860 
INT00870 
INT00880 
INT00890 
INT00900 
INTO0910 
INT00920 
INT00930 
INTO 09 40 
INT00950 
INT00960 
INT009  70 
INT00980 
INT00990 
INT01000 
INT01010 
INT01020 
INTO  1030 
INT01040 
INTO  1050 
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NRITE  =  (NI  +  D/2 

IMIN  =  (Nl-l)/2+l 

IF((Nl/2*2)  .  EQ.  Nl)  IMIN  =  Nl/2 

CALL  TRGRID  (Nl ,X0, YO ,NI , NRITE , 0. 5 , IMIN, RAD, 1 ,NXSS1) 

N1SAVE  =  Nl 
53    CONTINUE 

ALPHAA  =  0.0174533  *  ALPHAA 

ALPHAI  =  0.0174533  *  ALPHAI 

STAG   =  0. 0174533  *  STAG 
C 

IF  (INLET  .EQ.  0)  THEN 

ALPHA  =  ALPHAA 

ELSE 

ALPHA  =  ALPHAI 

END  IF 
C 

IF  (ISTAG  .NE.  0)  THEN 

CALL  STAGR( Nl , STAG , XO , YO , XSTGR , YSTGR) 

ELSE 

DO  55  I  =  1  ,  Nl 

XSTGR(I)  =  XO(I) 

YSTGR(I)  =  YO(I) 
55    CONTINUE 

END  IF 
C 

C  READ  DATA  FROM  VISCOUS  CAL. 
C 

I CYCLE  =  0 

ICYCLE  =  ICYCLE  +  1 


i 


60 
C 


70 
C 

999 


CALL  POTNL( N 1 , IRST , ALPHA , CHORD , XO , YO , XSTGR , YSTGR , X , Y , DLS , VCOM , 
+  DBPP) 

IF  (ICYCLE  .GT.  ICYTL)  THEN 
REWIND  3 

WRITE  (3)N1,(X0(I),Y0(I),DLS(I),VN(I),DBPP(I),I=1,N1) 
GOTO  1 
END  IF 

IF  (ISTAG  .NE.  0)  THEN 

DO  70  I  =  1  ,  Nl-1 

X(I)  =  0.5  *  (XO(I)+XO(I+l)) 

Y(I)  =  0.5  *  (YO(I)+YO(I+l)) 

CONTINUE 
END  IF 

CALL  CASBLP(Nl,XO,YO,X,Y,XS,YS, DLS, VCOM, DBPP, RN 
+  ,NBL,ITRI,XCTRI, TITLE) 

GO  TO  60 
CONTINUE 
STOP 
END 

SUBROUTINE  POTNL( Nl , IRST , ALPHA , CHORD , XO , YO , XSTGR , YSTGR , X , Y , DLS , 

+  VCOM, DBPP) 

COMMON/BLOW/VN(100) 


INT0106 

INT0107I 

INTO  1080 

INTO  1090 

INTO  1100 

INTO  11 10 

INTO  1120 

INTO  11 31 

INTO 1141 

INTO1150 

INT01160 

INTO 1171 

INTO  1180 

INTO1190 

INTO12O0 

INTO1210 

INT01220 

INTO1230 

INT01241 

INTO 1251 

INTO  1260 

INTO1270 

INTO1280 

INTO  1290 

INTO13O0 

INTO1310 

INTO1320 

INTO1330 

INTO1340 

INTO1350 

INTO  1360 

INTO 13 70 

INTO  1380 
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10 


15 


16 

17 
C 


20 


COMMON  /BLCO/  NX,NXT,NP ,NPT,NTR, IT, INVRS ,NS , IP 

COMMON/BLIN/TITLE(20) ,XC( 100) ,YC( 100) , ISG( 100) ,DELS( 100) , 
+  XCTR , XTR , ISTRP , ICYCLE , ICYTL,XCTRS( 2 ) ,TRFIND( 2 ) 

COMMON/CASCDE /  INLET , SP , SINGLE , ALPHAA , ALPHAI , STAG 
SIMPLE  SOURCE  POTENTIAL  CODE 

DIMENSION  AOFF(100,100),  EOFF( 100 , 100) ,  XP(IOO),  YP( 100) ,X( 100) , 
+     S(100),C( 100),  D(100),VTAN(3,100),VNOR(3,100),R(3,100), 
+     VCOM( 100) ,SIGCOM( 100) ,CP( 100) ,X0( 100) , Y0( 100) 
+    ,VNC(100),D1(100),D2(100),D3(100),SO(100),SC(100) 
+    ,XOFF( 100) ,YOFF( 100) ,T(3, 100) ,VTCOM( 100) ,VNCOM(  100) 
+     ,XS( 100) ,YS( 100) ,SOFF( 100) ,COFF( 100) ,XPOFF( 100) ,YPOFF( 100) 
+     ,Y(100),SIG(3,100),DLS(100),DLSC(100),A(100,100),B(100,100) 
+     ,XSTGR(100),YSTGR(100),VUT(3),VLT(3),VUN(3),VLN(3),DBPP(100) 
+     ,CPI( 100) ,XOS( 100) ,YOS( 100) ,DBPPC( 100) 

REAL  NUM1  ,  NUM2 

LOGICAL  OFF, SINGLE 

OFF  =  .FALSE. 

PI  =  3.  141592 

CM  =  0.0 

N  =  Nl  -  1 

IF  (ICYCLE  .EQ.  1)  THEN 

IF  (IRST  .EQ.  0)  THEN 

DO  10  1=1, Nl 

DLS(I)  =  0.  0 

VN  (I)  =  0.  0 

DBPP(I)=  0. 0 

CONTINUE 

ELSE 

DO  5  I  =  1  ,  Nl 

XOS(I)  =  XO(I) 

YOS(I)  =  YO(I) 

CONTINUE 

READ  (3)  NT,(XS(I),YS(I),DLSC(I),VNC(I),DBPPC(I),I=1,NT) 

XMIN  =  XS(1) 

DO  15  I  =  2  ,  NT 

IF  (XS(I)  .GT.  XMIN)  GOTO  15 

XMIN  =  XS(I) 

IMIN  =  I 

CONTINUE 

DO  17  I  =  1  ,  NT 

IF  (I  .LT.  IMIN)  GOTO  16 

XS(I)  =  XS(I)  -  XMIN 

GOTO  17 

XS(I)  =  XMIN  -  XS(I) 

CONTINUE 

XMIN  =  XOS(I) 

DO  20  I  =  2  ,  Nl 

IF  (XOS(I)  .GT.  XMIN)  GOTO  20 

XMIN  =  XOS(I) 

IMIN  =  I 

CONTINUE 

DO  22  I  =  1  ,  Nl 

IF  (I  .LT.  IMIN)  GOTO  21 

XOS(I)  =  XOS(I)  -  XMIN 
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21 
22 


GOTO  22 
XOS(I)  = 

CONTINUE 


XMIN  -  XOS(I) 


30 
C 


100 
C 


CALL  DIFF3(NT,XS,DLSC,D1,D2,D3,0) 
CALL  INTRP3(NT,XS,DLSC,D1,D2,D3,N1,X0S,DLS) 
CALL  AMEAN  ( 1 ,N1 ,XOS ,DLS , 1) 
CALL  DIFF3  (NT,XS , VNC ,D1 ,D2 ,D3 ,0) 
CALL  INTRP3(NT,XS,VNC,D1,D2,D3,N1,X0S,VN) 
CALL  AMEAN  ( 1 ,N1,X0S,VN, 1) 
CALL  DIFF3  (NT,XS ,DBPPC ,D1 ,D2 ,D3 ,0) 
CALL  INTRP3(NT,XS,DBPPC,D1,D2,D3,N1,X0S,DBPP) 
CALL  AMEAN  ( 1 ,N1,X0S,DBPP, 1) 
END  IF 
END  IF 

DO  30  1=1, Nl 
XP(I)  =  XSTGR(I) 
YP(I)  =  YSTGR(T) 
CONTINUE 
CALCULATE  GEOMETRIC  QUANTITIES 
DO  100  J=1,N 

VNC(J)  =0.5  *  (VN(J)  +  VN(J+1)) 
X(J)=  . 5*(XP(J)+XP(J+1)) 
Y(J)=  . 5*(YP(J)+YP(J+1)) 

D(J)  =  SQRT((XP(J+1)-XP(J))**2  +  (YP(J+1)-YP(J))**2) 
C(J)=  (XP(J+1)-XP(J))/D(J) 
S(J)=  (YP(J+1)-YP(J))/D(J) 


35 


NE.  0  .AND.  .NOT.  SINGLE)  THEN 

N 
D(  J) 


PI  *  SUM  /  SP 


C 
C 
102 


CONTINUE 

IF  (  INLET 
SUM  =  D(l) 
DO  35  J  =  2 
SUM  =  SUM  + 
CONTINUE 
Q  =  2.0 
ELSE 
Q  =  0.  0 
END  IF 

CALCULATE  NORMAL  AND  TANGENTIAL  MATRICES 
CONTINUE 

IF  (SINGLE)  THEN 
IF  (  .NOT.  OFF)  THEN 
DO  120  1=1, N 
DO  110  J=1,N 
IF  (J  . EQ.  I)  GO  TO  105 
XX=  (X(I)-X(J))*C(J)  +  (Y(I) 


YY=-(X(I)-X(J))*S(J)  +  (Y(I) 


Y(J))*S(J) 
Y(J))*C(J) 


105 


UU=  LOG( ( (XX+.  5*D( J) )**2+YY**2)/( (XX-.  5*D( J) )**2+YY**2)) 

VV=  2.  *ATAN2(YY*D( J) ,  XX**2+YY**2-(.  5*D( J)  )**2) 

SS=  S(I)*C(J)  -  C(I)*S(J) 

CC=  C(I)*C(J)  +  S(I)*S(J) 

A(I,J)=  -UU*SS  +  VV*CC 

B(I,J)=  UU-'-CC  +  W*SS 

GO  TO  110 

A(I,J)  =  6.2831853 
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B(I,J)  =  0. 0 
110   CONTINUE 
120   CONTINUE 

ELSE 

DO  140  1=1, N 

DO  130  J=1,N 

XX=  (XOFF(I)-X(J))*C(J)  +  (YOFF(I)-Y(J))*S(J) 

YY=-(XOFF(I)-X(J))*S(J)  +  (YOFF(I)-Y(J))*C(J) 

UU=  LOG( ( (XX+. 5*D( J))**2+YY**2)/((XX-.  5*D( J) )**2+YY**2)) 

VV=  2. *ATAN2( YY*D( J) ,  XX**2+YY**2 - ( .  5*D( J) ) **2 ) 

SS=  SOFF(I)*C(J)  -  COFF(I)*S(J) 

CC=  COFF(I)*C(J)  +  SOFF(I)*S(J) 


130 
140 


45 

40 
50 


•uu*ss 

UU*CC 


OFF)  THEN 


AOFF(I,J)= 

BOFF(I,J)= 

CONTINUE 

CONTINUE 

END  IF 

ELSE 

IF  (  .  NOT. 

DO  50  1=1, N 

DO  40  J=1,N 

IF  (J  .EQ.  I)  GO  TO  45 

XX=  (X(I)-X(J))*C(J)  + 

YY=-(X(I)-X(J)) 

DX1  =  PI  *(X(X) 

DY1  =  PI  *(YO) 

DX2  =  PI  *(X(I) 

DY2  =  PI  *(Y(I) 

risq  =  (coswrn-: 

(COSH(DX2)) 

L0G(R1SQ/R2SQ) 


vv*cc 
vv*ss 


R2SO  = 

uu 

NUM1  = 


'S(J)  + 

XP(J)) 

YP(J)) 

XP(J+1)) 

YP(J+1)) 

))**2 


(Y(I) 
(Y(I) 

/  SP 
/  SP 

/  SP 
/  SP 


Y(J))*S(J) 
Y(J))*C(J) 


+ 


DNUM1= 


+ 


NUM2  = 


DNUM2= 


DX1 
DY1 
DX1 
DY1 
DX2 
DY2 
DX2 
DY2 
2.  0 


COSH(DXl) 
SINH(DXl) 
SINH(DXl) 
COSH(DXl) 
COSH(DX2) 
SINH(DX2) 
SINH(DX2) 
COSK(DX2) 


(C0S(DY1))**2 
(C0S(DY2))**2 

r  SIN(DYl)  - 

•  COS(DYl) 

■  COS(DYl)  + 

■  SIN(DYl) 

'  SIN(DY2)  - 
'  C0S(DY2) 

•  COS(DY2)  + 

•  SIN(DY2) 


ATAN2(NUM2,DNUM2)  ■ 
XX**2+YY** 


EXV  = 

VV=  2 . *ATAN 2 ( YY*D ( J ) 

VV  =  VV  +  EXV 

SS=  S(I)*C(J)  -  C(I)*S(J) 

CC=  C(I)*C(J)  +  S(I)*S(J) 

A(I,J)=  -UU*SS  +  VV*CC 

B(I,J)=  UU*CC  +  VV*SS 

GO  TO  40 

A(I,J)  =  6. 2831853 

B(I,J)  =0.0 

CONTINUE 

CONTINUE 

ELSE 

DO  70  1=1, N 

DO  60  J=1,N 

XX=  (XOFF(I)-X(J))*C(J) 

YY=-(XGFF(I)-X(J))*S(J) 


2.0  *  ATAN2(NUM1,DNUM1) 
2-(.5*D(J))**2) 


+  (YOFF(I)-Y(J))*S(J) 
+  (YOFF(I)-Y(J))-C(J) 
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YP(J)) 


DX1  =  PI  *(XOFF(I)-XP(J)) 
DY1  =  PI  *(YOFF(I) 
DX2  =  PI  *(XOFF(I) 
DY2  =  PI  *(YOFF(I) 
R1SQ  =  (COSH(DXl)) 
R2SQ  =  (COSH(DX2))**2 
=  L0G(R1SQ/R2SQ) 


/ 
/ 


XP(J+1)) 
YP(J+1)) 


SP 

SP 

/  SP 
/  SP 


1**2 


(COS(DYl))**2 
(COS(DY2))**2 


UU 
NUM1  = 

DNUM1= 

h 

NUM2  = 
h 
DNUM2= 


+ 


DX1 
DY1 
DX1 
DY1 
DX2 
DY2 
DX2 
DY2 
2.0 


* 
* 


COSH(DXl) 
SINH(DXl) 
SINH(DXl) 
COSH(DXl) 
COSH(DX2) 
SINH(DX2) 
SINH(DX2) 
C0SH(DX2) 
ATAN2(NUM2,DNUM2) 
XX**2+YY* 


SIN(DYl) 
COS(DYl) 
COS(DYl) 
SIN(DYl) 
SIN(DY2) 
COS(DY2) 
COS(DY2) 
SIN(DY2) 


60 
70 


C 

c 


145 


EXV  =  2.0  *  ATAN2(NUM2,DNUM2)  -  2.  0  ' 
VV=  2 . * ATAN2 ( YY*D( J ) ,  XX**2+YY**2 - ( .  5 

VV  =  VV  +  EXV 

SS=  SOFF(I)*C(J)  -  COFF(I)*S(J) 
CC=  COFF(I)*C(J)  +  SOFF(I)*S(J) 
AOFF(I,J)=  -UU*SS  +  VV*CC 
BOFF(I,J)=  UU*CC  +  VV*SS 
CONTINUE 
CONTINUE 
END  IF 
END  IF 
NORMAL  AND  TANGENTIAL  COMPONENTS  OF 

DO  160  1=1, N 

SUMR=  0. 

SUMT=  0. 

IF  (   .NOT.  OFF)  THEN 

R(1,I)=  S(I)+VNC(I)/COS(ALPHA) 

T(1,I)=  C(I) 

R(2,I)=  -C(I) 

T(2,I)=  S(I) 

DO  145  J=1,N 

SUMR  =  SUMR  +  B(I,J) 

SUMT  =  SUMT  +  A(I,J) 

CONTINUE 

ELSE 


*  ATAN2(NUM1,DNUM1) 
'D(J)),V*2) 


FUNDAMENTAL  SOLUTIONS 


R(1,I) 

T(1,I) 
R(2,I) 
T(2,I) 
DO  150 

SUMR= 


=  SOFF(I) 
=  COFF(I) 
=-COFF(I) 
=  SOFF(I) 
J=1,N 
SUMR  + 


150 


160 
C 


SUMT=  SUMT 
CONTINUE 
END  IF 

R(3,I)=  SUMR 
T(3,I)=  SUMT 
CONTINUE 

IF  (  OFF  )  GO 

DECOMPOSITION 


BOFF(I,J) 
AOFF(I,J) 


TO 
OF 


275 

MATRIX 
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DO  230  1=1, N-l 

DO  220  K=I+1,N 

A(K,I)=  A(K,I)/A(I,I) 

DO  210  J=I+1,N 

A(K,J)=  A(K,J)-  A(K,I)*A(I,J) 

210 

CONTINUE 

220 

CONTINUE 

230 

CONTINUE 

C   OPERATE  ON  FUNDAMENTAL  RIGHT  SIDES 

DO  270  K=l,3 

DO  260  J=1,N-1 

DO  250  I=J+1,N 

R(K,I)=  R(K,I)  -  A(I,J)*R(K,J) 

250 

CONTINUE 

260 

CONTINUE 

270 

CONTINUE 

C   BACK  SOLUTION 

DO  300  K=l,3 

DO  290  I=N,1,-1 

SUM=  0. 

DO  280  J=N, 1+1,-1 

SUM=  SUM  +  A(I,J)*SIG(K,J) 

280 

CONTINUE 

SIG(K,I)=  (R(K,I)-SUM)/A(I,I) 

290 

CONTINUE 

300 

CONTINUE 

OFF  =  .TRUE. 

C 

c 

ADD  DIS -PLACE  VERTICALLY  TO  THE 

c 

DISPLACEMENT  SURFACE 

c 

WITH  LOWER  TRIANGULAR 


BODY  TO  GENERATE 


1)*D(I))/(D(I-1)+D(I)) 
1)*D(I))/(D(I-1)+D(I)) 


DO  305  I  =  2  ,  N 

SOFF(I)  =  (S(I)*D(I-1)+S(I' 

COFF(I)  =  (C(I)*D(I-1)+C(I- 

305  CONTINUE 

SOFF(l)  =  2. 0*S(1)  -  SOFF(2) 

S0FF(N1)=  2. 0*S(N)  -  SOFF(N) 

COFF(l)  =  2.0-C(l)  -  COFF(2) 

C0FF(N1)=  2.0*C(N)  -  COFF(N) 

DO  306  I  =  1  ,  Nl 

XPOFF(I)  =  XP(I)  -  SOFF(I)  *  DLS(I) 

YPOFF(I)  =  YP(1)  +  COFF(I)  *  DLS(I) 

306  CONTINUE 

DO  307  I  =  1  ,  N 

XOFF(I)  =  0.5  *  (XPOFF(I)  +  XPOFF(I+l)) 
YOFF(I)  =0.5  *  (YPOFF(I)  +  YPOFF(I+l)) 
DOFF  =  SQRT((T'POFF(I+l)-XPOFF(I))**2  + 
+  (IPOFF(I+l)-YPOFF(I))**2  ) 

COFF(I)  =  (XPOFF(I+l)-XPOFF(I))/DOFF 
SOFF(I)  =  (YPO^(I  +  l)-YPOFF(I))/DOFF 
CONTINUE 
GO  TO  102 
CALCULATION  OF  SURFACE  VELOCITIES  FOR  THE  FUNDAMENTAL  SOLUTIONS 


307 

C 
C 
275 


DO 

DO 


330 
320 


K=l 

1  =  1 
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310 


320 
C 


SUMT=  T(K,I) 

SUMN=-R(K,I) 

DO  310  J=1,N 

SUMT=  SUMT  +  BOFF(I,J)*SIG(K,J) 

SUMN=  SUMN  +  AOFF(I,J)*SIG(K,J) 

CONTINUE 

VTAN(K,I)=  SUMT 

VNOR(K,I)=  SUMN 

CONTINUE 


D0FF1  =  SQRT( (XPOFF( 2) -XPOFF( 1) )**2+( YPOFF( 2) -YPOFF( 1) )**2) 
DOFF2  =  SQRT( (XPOFF( 3) -XPOFF( 2) )**2+( YPOFF( 3) -YPOFF( 2) )**2) 

DOFFN  =  SQRT((XPOFF(N+l)-XPOFF(N))**2+(YPOFF(N+l)-YPOFF(N)) 
+  **2) 

D0FFN1=  SQRT((XPOFF(N)-XPOFF(N-l))**2+(YPOFF(N)-YPOFF(N-l)) 
+  **2) 

VUT(K)  =  VTAN(K,N)  +  DOFFN  *  ( VTAN(K,N) -VTAN(K,N-1) )/ 
+  (DOFFN+DOFFNT) 

VUN(K)  =  VNOR(K,N)  +  DOFFN  *  ( VNOR(K,N) -VNOR(K,N-l) )/ 
+  (DOFFN+DOFFN1) 

VLT(K)  =  VTAN(K,1)  +  DOFF1  *  ( VTAN(K, 1) -VTAN(K, 2) )/ 
+  (DOFF1+DOFF2) 

VLN(K)  =  VNOR(K,l)  +  DOFF1  *  ( VNOR(K, 1) -VNOR(K, 2) )/ 
+  (DOFF1+DOFF2) 

330    CONTINUE 
C   OUTPUT  FUNDAMENTAL  SOLUTIONS 

IF  (ICYCLE  .EQ.  1  .OR.  ICYCLE  .  GE.  ICYTL-1  .OR.  IP  .  GE.  0) 
+  WRITE(6,335)  TITLE 
335   F0RMAT(1H1,///20A4//) 
C     DO  360  K=l,3 
C     WRITE  (6,340)  K 

C340   F0RMAT(////,1H  ,' FUNDAMENTAL  SOLUTION  NUMBER  ',12////) 
C      WRITE(6,345) 

345   F0RMAT(3X,'I' ,8X, ' X' , 11X, ' Y1 , 10X, ' VT1 , 10X, ' VN' ,8X, ' SIGf  ///) 
C      WRITE(6,375)  1,   XP(1)  ,  YP(1),  VLT(K) ,VLN(K) 
C      DO  350  1=1, N 

C      WRITE(6,375)  I  ,  X(I),  Y(I),  VTAN(K, I) , VNOR(K, I) ,SIG(K, I) 
C350    CONTINUE 

C      WRITE(6,375)  Nl  ,  XP(N1) , YP(N1) ,VUT(K) ,VUN(K) 
C360    CONTINUE 

C   COMBINED  FLOW  AT  ANGLE  OF  ATTACK 
C 

IF  (  INLET  .NE.  0)  THEN 

YYY  =  ((VUT(3)+VLT(3))'VTAN(ALPHAI)+(VUT(1)+VLT(1))^Q)/ 
+      ((VTJT(3)+VLT(3))-(VUT(2)+VLT(2))*Q) 

XXX  =  -((VUT(1)+VLT(1))+(VUT(2)+VLT(2))*TAN(ALPHAI))/ 
+       ((VUT(3)+VLT(3))-(VUT(2)+VLT(2))*Q) 
ALPHA  =  ACOS(l.  0/SQRT(l.  O+YYY—2)) 
COSAL  =  COS( ALPHA) 
SINAL  =  SIN (ALPHA) 
W  =  XXX/SQRT(1.0+YYY**2) 
ELSE 

COSAL  =  COS( ALPHA) 
SINAL  =  SIN( ALPHA) 

W=-((VLT(l)+VUT(l))*COSAL+(VLT(2)+VUT(2))*SINAL)/ 
+       (VLT(3)+VUT(3)) 


INT04410 
INT04420 
INT04430 
INT04440 
INT04450 
INT04460 
INT04470 
INT04480 
INT04490 
INT04500 
INT04510 
INT04520 
INT04530 
INT04540 
INT04550 
INT04560 
INT04570 
INT04580 
INT04590 
INT04600 
INT04610 
INT04620 
INT04630 
INT04640 
INT04650 
INT04660 
INT04670 
INT04680 
INT04690 
INT04700 
INT04710 
INT04720 
INT04730 
INT04740 
INT04750 
INT04760 
INT04770 
INT04780 
INT04790 
INT04800 
INT04810 
INT04820 
INT04830 
INT04840 
INT04850 
INT04860 
INT04870 
INT04880 
INT04890 
INT04900 
INT04910 
INT04920 
INT04930 
INT04940 
INT04950 
INT04960 


72 


c 
c 
c 


END  IF 
C   FORCE  COEFFICIENT  CALCULATION 

SUM1=  0. 

SUMX=  0. 

SUMY=  0. 

DO  390  1=1, N 

SUM1=  SUM1+  D(I) 

SUMX=  SUMX-  VCOM( I)**2*S( I  )*D( I ) 

SUMY=  SUMY+  VCOM( I)**2*C( I)*D( I ) 
390   CONTINUE 
C  FIND  MAN.  CHORD  LENGTH 

XOMIN  =  XO(l) 

DO  395  I  =  2  ,  Nl 

IF  (  XO(I)  . GT.  XOMIN)  GOTO  395 

XOMIN  =  XO(I) 
395    CONTINUE 

CHORD  =  XO(Nl)  -  XOMIN 

CL1=  SUMl-->25.  13274*W/CHORD 

CL2=  ( SUMY*COSAL-SUMX*SINAL) /CHORD 

CD  =  (SUMX*COSAL+SUMY*SINAL) /CHORD 

CALCULATING  PARAMETERS  FOR  INLET  VELOCITY  AS  MODULUS  OF  NOMORIZED 

IF  (.NOT.  SINGLE)  THEN 

NUM1  =  SIN(ALPHA)+CLl*CH0RD/(4.  0*SP) 

ALPHID  =  ATAN2(NUM1,C0S( ALPHA)) 

NUM1  =  SIN( ALPHA) -CLl*CH0RD/(4.  0*SP) 

ALPHED  =  ATAN2(NUM1, COS (ALPHA)) 

NUM1  =  CLl*CHORD/(2.  0*SP)*COS(ALPHA) 

DNUM1=  1.  0-(CLl*CH0RD/(4.  0*SP))**2 

DALPKA  =  ATAN2(NUM1,DNUM1) 

UOUI  =  (TAN(ALPHID)-TAN(ALPHED))*(2.  0*SP/CHORD*COS(  ALPHID)  )/CLl 

CLI  =  CL1  *  UOUI**2 

UIOU  =  1. 0/UOUI 

VEXIT  =  COS(ALPHA)/COS(ALPHED) 

ELSE 

ALPHID  =  ALPHA 

ALPHED  =  ALFHA 

DALPHA  =  0.0 

UOUI  =  1.0 

UIOU  =1.0 

CLI   =  CLI 

VEXIT  =1.0 

END  IF 

FAC  =  180. 0/PI 

ALPHID  =  ALPHID 

ALPHED  =  ALPHED 

DALPHA  =  DALPHA 

ALPHAD  =  ALPHA 

IF  (ICYCLE  .  EQ. 


370 


*  FAC 

*  FAC 

*  FAC 

*  FAC 

1  .OR.  ICYCLE  .  GE.  ICYTL-1  .OR.  IP  .  GE.  0)  THEN 
WRITE(6,370)  ALPHAD  ,ALPHID, ALPHED, DALPHA, UIOU, VEXIT 
FORMAT  (////, IK  ,  'COMBINED  FLOW  AT  AVERAGE  ANGLE  OF  ATTACK  =  ' 
+  F8.3,  4X, 'DEGREES' ,  / , 1H  . 17X, ' INLET  ANGLE  OF  ', 

+         'ATTACK   =  ' ,F8. 3, 4X, 'DEGREES' ,/,lH  , 
+  17X, 'EXIT  ANGLE  =   * ,F8. 3 ,4X, ' DEGREES' ,/, 1H  ,17X, 


TURNNING  ANGLE  =  ' ,F8. 3 ,4X, ' DEGREES ',/, 1H  ,17X, 


INT04970 
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+  'INLET  VEL  =  ' ,F10. 6,3X,'EXIT  VEL  =  f,F10. 6,//) 

WRITE(6,365) 
365    F0RMAT(3X, ' i' ,8X, 'XO' ,10X, 'YO' ,10X, 'X! ,11X, *Y' ,10X, 'VT' , 
+  10X,'VN'  ,11X,V  ,10X,'CP'  ,9X,'CPI'///) 

END  IF 

DO  380  1=1, N 

VTCOMC I )=  VTAN( 1 , I )*COSAL+VTAN( 2 , I ) *SINAL+W*VTAN( 3,1) 

VNCOM( I )=  VNOR( 1 , I )*COSAL+VNOR( 2 , I )*SINAL+W*VNOR( 3,1) 

VCOM(I)=  SQRT( VTCOM( I )**2  +  VNCOM(I)**2) 

IF  (VTCOM( I)  . LT.  0.0)  VCOM(I)  =  -VCOMCI) 

CP(I)  =  1.0  -  VCOMCI)  **  2 

CPI(I)=  1.0  -  (VCOM(I)*UOUI)**2 

SIGCOM(I)  =  SIG(l,I)*COSAL+SIG(2,I)*SINAL+W*SIG(3,I) 

XP(I)  =0.5  *  CXO(I)+XO(I  +  l)) 

YP(I)  =  0.5  *  (YO(I)+YO(I+l)) 
380   CONTINUE 

IF  (ICYCLE  .  EQ.  1  .OR.  ICYCLE  .  GE.  ICYTL-1  .OR.  IP  .  GE.  0)  THEN 

WRITE  (1,374)   C  XO(I),YO(I),XP(I),YP(I),CP(I),CPI(I)  ,1=1, N) 

WRITE  (6,375)   (  I,  XO(I),  YO(I),  XP(I),  YP(I),  VTCOMC I), 
+  VNCOM(I),VCOM(I),CP(I),  CPI(l)  ,I=1,N) 

FORMAT(6F10.4) 

FORMATC1X,  13,  9F12. 4) 

WRITE  (2)  I,XO(I),YO(I),XSTGR(I),YSTGR(I),DLS(I),X(I),Y(I),VCOM( 

WRITE(6,385)  N+l ,XO(N+l) , YO(N+l) 

F0RMAT(1X,I3,2F12. 4) 

WRITE(6,400)  CHORD,  CL1  ,CLI 

FORMAT(///3X,' CHORD  =  ' ,F10. 5 ,4X, ' CLC AVG)  =  ',F10.5,4X, 
+  'CL( INLET)  =   ,F10.5) 

END  IF 

F0RMAT(/3X,1HI,6X,2HS0,10X,2HSC,9X,3H  VN,9X,3HVNC,9X,3HDLS ,8X, 
+  4HDLSC) 

FORMATCI5,6E12.4) 

RETURN 

END 


374 
375 
C 

385 

400 

420 
430 


C 
C 
C 
C 

C 

c 
c 
c 
c 
c 
c 
c 

c 
c 
c 


DATA  SET  KCBCAMEAN  AT  LEVEL  001  AS  OF  08/24/84 
DATA  SET  KCBCAMEAN  AT  LEVEL  003  AS  OF  04/05/84 
SUBROUTINE  AMEAN(NS,ND,X,Y, IT) 

SMOOTH  DATA  USING  3-PTS  WEIGHTING  FORMULA 


NS 
ND 
X,  Y 


IT 


STARTING  PINT  OF  THE  DATA  TO  BE  SMOOTHED 
END  PINT  OF  THE  DATA  TO  BE  SMOOTHED 
INDEPENDENT  +  DEPENDENT  VARAIBLES  OF  THE  DATA 
TO  BE  SMOOTHED 
CYCLES  OF  DATA  SMOOTHING 


DIMENSION  X(101),Y(101) 


NM    =  ND  -NS 

IF(NM  .LT.  2  .OR.  IT  .  LT.  1)  RETURN 

NDM1  =  ND  -  1 
NSP1  =  NS  +  1 
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DO  20  K=1,IT 

INT06090 

DL1    =  X(NSP1)  -  X(NS) 

INT06100 

Yl    =  Y(NS) 

INT06110 

DO  10  I=NSP1,NDM1 

INT06120 

DL2   =  X(I  +  1)  -X(I) 

INT06130 

Y2    =  Y(I) 

INT06140 

YM    =  (DL2  ■'  Yl  +  DL1  *  Y(I+1))/(DL] 

.  + 

DL2) 

INT06150 

Y(I)   =  0. 5  *  (Y2  +  YM) 

INT06160 

DL1   =  DL2 

INT06170 

Yl    =  Y2 

INT06180 

10 

CONTINUE 

INT06190 

20 

CONTINUE 

INT06200 

C 

RETURN 
END 

INT06210 
INT06220 
INT06230 

C 

DATA  SET  KCBCBLGRID  AT  LEVEL 

001 

AS 

OF 

08/24/84 

INT06240 

C 

DATA  SET  KCBCBLGRID  AT  LEVEL 

001 

AS 

OF 

08/24/84 

INT06250 

C 

DATA  SET  KCBCBLGRID  AT  LEVEL 

004 

AS 

OF 

04/05/84 

INT06260 

SUBROUTINE  BLGRID(N,X,T,D1) 

INT06270 

c 

INT06280 

c 

GENERATE  B.  L.  X-V'ISE  GRID  USING  MODI 

DISTRIBUTION 

INT06290 

c 

DIMENSION  X(101),T(101),D1(101) 

INT06300 
INT06310 

DATA  CRAD/57. 2957795/,  BPI/3. 14159265/ 

r 

INT06320 

c 

INT06330 

c 
c 

INT06340 
INT06350 

NN     =2  *  N  -  1 

INT06360 

EN     =  FLOAT((NN-l)/2) 

INT06370 

THO    =  10. /CRAD 

INT06380 

CTOl   =  1.  +  COS(THO) 

INT06390 

DTH    =  (BPI  -  THO)  /  EN 

INT06400 

FI     =  FLOAT(N  -  2) 

INT06410 

DO  10  I=N,NN 

INT06420 

FI     =  1.  0  +  FI 

INT06430 

II     =  I  -  N  +  1 

INTO 6 440 

XII    =  THO  +  FI  *  DTH 

INT06450 

X(II)   =  (1.0  +  C0S(XII))/CT01 

INT06460 

10 

CONTINUE 

XI     =  X(l) 

XN     =  X(N) 

CH     =  XN  -XI 

FN1    =  FLOAT(N-l) 

N10     =  N/10 

DO  20  1=1, N 

T(I)   =  FL0AT(I-1)/FN1 

X(I)    =  (X(I)-X1)/CH 

INT06470 
INTO 6480 
INT06490 
INT06500 
INT06510 
INT06520 
INT06530 
INT06540 
INT06550 

20 

CONTINUE 

INT06560 

c 

CALL  SMFIT(N10,N,T,X,D1,N10) 

INT06570 

c 

IF(X(2).LT.  0.  35  *  X(3))  X(2)  =  0.: 

35  * 

X(3) 

INT06580 

c 

CALL  SMFIT(1,N,T,X,D1,2) 
CALL  AMEAN(1,N,T,X,N10) 

INT06590 
INT06600 

c 

RETURN 
END 

INT06610 
INT06620 
INT06630 

c 

DATA  SET  KCBCBL2D   AT  LEVEL 

001 

AS 

OF 

08, 

Z24/84 

INT06640 
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c 
c 

c 
c 
c 


DATA  SET  KCBCBL2D  AT  LEVEL  001  AS  OF  08/24/84 
DATA  SET  KCBCBL2D  AT  LEVEL  012  AS  OF  04/06/84 
SUBROUTINE  BL2D  (  ITR , ISWPT, SURFID) 

PROGRAM  CALCULATES  VISCOUS/INVISCID  INTERACTION  USING  HILBERT 

INTEGRAL. 


C 
C 

c 
c 
c 


c 
c 
c 

c 
c 
c 


c 

10 


20 


30 


COMMON  /BLCO/ 
COMMON  /BLC1/ 
COMMON  /BLC2/ 
COMMON  /BLC7/ 
C0MM0N/EDDY1/ 


NX,NXT,NP,NPT,NTR,IT,INVRS,NS,IP 
F(101,2),U(10i,2),V(101,2),W(101,2),B(101,2) 
DELF( 101) ,DELU( 101) ,DELV( 101) ,DELW( 101) 
C(100,100),D(100)JDB(100),DBP(100),UE0(100)JGI 
RL,RX,SQRX,RXNTR,GMTR,GMTRS( 100) ,ALFAS( 100) , 
+  FFS(IOO) ,RTS(100),IEDY,NXSPT 

COMMON  /SMRY/  W( 100) , ITP( 100) , ISL( 100) ,DLS( 100) ,CF( 100) , 
+  THT(100),NPSTR(100) 

COMMON  /GTY  /  X( 101) ,UE( 100) ,P1( 100) ,P2( 100) ,CEL,CELH 

COMMON  /BONV/  ITMAX,EPSL,EPST,CONV 

COMMON  /SAVE/  FS( 101) ,US( 101) , VS( 101) ,WS( 101) ,BS( 101) 

COMMON  /BLIN/  TITLE( 20) ,XC( 100) , YC( 100) , ISG( 100) ,DELS( 100) , 
+  XCTR,XTR,ISTRP,ICYCLE,ICYTL,XCTRS(2),TRFIND(2) 

COMMON  /1SURF/  ISF 

COMMON/PLOT /v": ,2) ,NXVP( 20 , 2) , ICC 

DIMENSION  SURFID(4) 


GENERATE  B.  L.  GRIDS  +  SET  INITIAL  CONDITIONS 

DO  5  I  =  1  ,  NXT 

ALFAS(I)  =  0.  0 

FFS(I)    =  1.  0 

RTS(I)    =  1.  0 

CONTINUE 

CALL  INPUT( ITR, ISWPT, SURFID) 

CALCULATE  HILBERT  COEFFS. ,  C(I,J) 

CALL  CALCIJ(NXT,0) 

LOOP  OF  CALCULATIONS 

NSS   =  NS 
NXSPT  =  NXT  +  1 

IF  (  ICYCLE  .  EQ.  1  )  NS  =  NXT  +  1 
NX    =  NX  +  1 

CEL   =0.5  *  (X(NX)  +  X(NX-l))  /(X(NX)  -X(NX-l)) 
CELH  =0.5  *  CEL 
IT    =0 

RX    =  UE(NX)*X(NX)*RL 
SQRX  =  SQRT(RX) 
IT    =  IT  +  1 
IF(IT  . LE.  ITMAX)  GO  TO  40 
NXM1   =  NX-1 

CALL  HEADER(  TITLE , SURFID, ISTRP  ) 

WRITE(6,  170  )  (M,X(M),CF(M),DLS(M),UE(M),P2   (M),THT(M), 
+  D(M),ALFAS(M),   ITP(M) ,NPSTR(M) ,M=1 ,NXM1) 

WRITE(6,  160  )  NX 
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c 
c 
c 
c 


60 


70 

80 
C 
C 
C 


STOP 
40    CONTINUE 

IFCNX  .GT.  NTR)  CALL  EDDY 
CALL  COEFTR 
CALL  SOLV3 
IF(V(1,2).GT.  0.  0)  GOTO  60 

EXTRAPOLATE  CALCULATED  D  FOR  TURBULENT  SEPARATION  OR  LAMINAR 
SEPARATION  FOR  LAMINAR  FLOW  CALCULATION  ONLY 

CALL  EXTRAP(NX,NXT,X,D) 
NXM1  =  NX  -  1 

CALL  HEADER(  TITLE , SURFID , ISTRP  ) 

WRITE(6,  170  )  (M,X(M),CF(M),DLS(M),UE(M),P2(M),THT(M), 
+  D(M),ALFAS(M),ITP(M),NPSTR(M),M=1,NXM1) 

WRITE(6,180) 

WRITE(6,190)  (M,X(M),D(M),M=NX,NXT) 
GOTO  130 

IF(NX  .GT.  NTR)  GO  TO  70 
IF(ABS(DELV(lj)  . GT.  EPSL)  GO  TO  30 
GO  TO  SO 
CONTINUE 

IF(ABS(DELV(1)/V(1,2))  .  GT.  EPST)  GO  TO  30 
CONTINUE 

CHECK  FOR  GROWTH 

IF(NP  . GE.  NP7)  GO  TO  90 

IF(ABS(V(NP,2;;  .  LT.  0.0005  .AND.  ABS(  1.  0-U(NP-2 ,2)) 
+  .LT.  0.  0035)  GOTO  90 

CALL  FILLUP(l) 
IT  =  1 
GO  to  30 
CONTINUE 

CALL  FILLUP(2) 

CALL  OUTPUT(l) 

IF(ITR.EQ.  0  .OR.  NX.  GE.  NTR)  GOTO  100 

IF(NX.  LT.  3  .OR.  ITR.NE.  3)  GO  TO  100 

CALL  TRNS(ICODE) 

IF(ICODE.EQ.  1)  GOTO  20 

IF(NX  .NE.  NSS)  GOTO  120 

STORE  PROFILES  AT  THE  STATION  NS  FOR  INVERSE  B.  L. 
CALCULATION 

DO  110  J  =  1  ,  NPT 

FS(J)  =  F(J,2) 

US(J)  =  U(J,2) 

VS(J)  =  V(J,2) 

WS(J)  =  W(J,2) 

BS(J)  =  B(J,2) 
110   CONTINUE 

120    IF  (  NX  . LT.  NSS  )  GOTO  10 
C     IF  (ICYCLE  .NE.  1)  GO  TO  130 

IF  (  NX  .GE.  NS  )   GOTO  130 


90 
C 


100 

c 
c 
c 
c 
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INT07750 
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130 
140 


150 
155 
C 
C 

C 
C 
C 


158 

C 

C 

c 

160 


180 
190 
C 

n 


c 
c 
c 


c 
c 
c 


IF  (NX  .  LT 
CALL  KEADER( 
WRITE(6,  170 


+ 


NXT)  GO  TO  10 
TITLE,SURFID,ISTRP  ) 

)  (M,X(M),CF(M),DLS(M),UE(M),P2   (M),THT(M), 
D(M)5  ALFAS(M),  ITP(M) ,NPSTR(M) ,M=1 ,NXT) 
,  NXT 


NPT 


DO  140  I  =  1  , 
DB(I)  =  D(I) 
NS  =  NSS 
NX  =  NS 
NP  =  NPSTR(NX) 
DO  150  J  =  1  , 
F(J,2)  =  FS(J) 
U(J,2)  =  US(J) 
V(J,2)  =  VS(J) 
W(J,2)  =  WS(J) 
B(J,2)  =  BS(J) 
CONTINUE 
INVRS  =  NS  +  1 


CALCULATION  SHIFTS  TO  USING  PHYSICAL  COORDINATES 
CALL  MATN2(ITR,ISWPT,SURFID) 

PASS  DELTA-STAR  BACK  TO  MAIN  PROG. 

DO  158  I  =  1,NXT 
DELS(I)  =  DLS(I) 
CONTINUE 
RETURN 


FORMAT(lH0,'  **  ITERATIONS  EXCEEDED  ITMAX  AT  NX  =  ',15/ 


1H  , '  **  CALCULATIONS  STOP.  **  ' ) 
**  SUMMARY  OF  STANDARD  B.  L.  SOLUTIONS.  **'/ 
1H0,4X,2HNX,7X,1HX,12X,2HCF,11X,3HDLS,12X,2HUE, 
12X,2HP2,11X,3HTHT,13X,1HD,10X,4HALFAS6X,2HIT,2XS 
(1H  ,3X,I3,F10.5,2X,7E14.5,2I4)) 
FORMAT(1HO,34H  FLOW  SEPARATES.  D  IS  EXTRAPOLATED/ 
1H0,3X,3H  NX,7X,1HX,13X,1HD/) 
,3XI3,F10.5,2X,E14.  5) 


170   F0RMAT(1H0 

+ 
+ 
+ 

+ 


FORMAT( 1H 

END 

DATA 
DATA 
DATA 

SUBROUTINE 


SET  KCBCCALCIJ  AT  LEVEL  001  AS  OF  08/24/84 

SET  KCBCCALCIJ  AT  LEVEL  001  AS  OF  08/24/84 

:L7  KCBCCALCIJ  AT  LEVEL  005  AS  OF  04/05/84 
CALCIJ  (  IL,  LO) 


CALCULATE  HILBERT  INTEGRAL  COEFFS 

COMMON  /BLC7/  C( 100 , 100) ,D( 100) ,DB( 100) ,DBP( 100) ,UE0( 100) ,GI 
COMMON/EDDY1/  RL,RX, SQRX,RXNTR,GMTR,GMTRS( 100) 
+  ,ALFAS( 100) ,FFS( 100) ,RTS( 100) , IEDY,NXSPT 

COMMON  /GTY  /  X( 101) ,UE( 100) ,P1( 100) ,P2( 100) ,CEL,CELH 
DIMENSION    E(2) 


PI 


=  3. 14159265 


INTO 7 77 
INTO 7 78 
INTO 7 79 
INTO 7 8^ 
INTO 7 8 
INTO 7 8 
INT078 
INT078 
INT078 
INT0786 
INT0787 
INTO 7 88 
INTO 7 89 
INT0790 
INT079 
INT079 
INT0793 
INTO 7 94 
INT0795 
INT0796 
1NT0797 
INT0798 
INT0799' 
INT08001 
INTO  8  01' 
INTO  802i 
INT0803I 
INT0804I 
INT0805I 
INT0806I 
INT0807I 
INT0808I 
INT0809I 
INT0810I 
2HNP/INT0811! 
INTO  8 121 
INTO  8 131 
INT0814I 
INT0815I 
INT0816I 
INT0817 
INTO  8 18i 
INT0819' 
INTO  8  20' 
INT0821' 
INT08221 
INT0823- 
INT0824 
INT0825 
INT0826 
INT0827' 
INT0828 
INT0829 
-  -  INT0830 
INT0831 
INT0832 


78 


c 
c 

c 


c 
c 

c 

30 


c 
c 
c 

40 

c 

50 


60 


65 

C 

c 


c 
c 
c 

c 
c 
c 
c 


PI      =  PI*SQRT(RL) 

IL1     =  IL   -  1 

DO  65  I  =  2,  IL1 

E  (1)   =  0. 

L      =  LO  +  I 

DO  60  J  =  2,  IL 

Jl      =  J  -  1 

K       =  J  +  LO 

DX1     =  X(L)  -  X(K) 

DX2     =  X(K)  -  X(K-l) 

DX3     =  X(L)  -  X(K-l) 

IF  (  J  .  EQ.  I  )  GO  TO  30 

IF  (  J  .EQ.  (1  +  1)  )  GO  TO  40 

J  .NE.  I  OR  1  +  1 

E  (2)   =  (  1.0/DX2  )  *  ALOG(  ABS(  DX3  /  DX1  )  ) 

GO  TO  50 

J  .EQ.  I 

Rl      =  (  X(K+1)-X(K)  )  /  (X(K+1)-X(K-1)  ) 

E  (2)   =  (  Rl  *  ALOG(  ABS(  DX3  /  (  X(L) -X(K+1) )  )  )  +  2.  0 

GO  TO  50 

J  .EQ.  1  +  1 

Rl      =  (  X(K-l)-X(K-2)  )  /  (  X(K)-X(K-2)  ) 

E  (2)   =  (  Rl  *  ALOG(  ABS(  (X(L)-X(K-2) )  /  DX1  )  )  -  2.  0 

CONTINUE 

C  (J1,I)=  (  E(l)  -  E(2)  )  /  PI 

E(l)    =  E  (2) 

CONTINUE 

E  (2)    =  0. 

Jl      =  IL 

K       =  K  +  1 

C  (J1,I)=  E(l)  /  PI 

CONTINUE 


RETURN 

END 

DATA  SET  KCBCCOEF 
DATA  SET  KCBCCOEF 
DATA  SET  KCBCCOEF 


AT  LEVEL  001  AS  OF  08/24/84 
AT  LEVEL  001  AS  OF  08/24/84 
AT  LEVEL  007  AS  OF  04/05/84 


SUBROUTINE  COEF( GAMMA 1 ,GAMMA2) 

CALCULATE  COEFFS  OF  B.  L.  FINITE-DIFFERENCE  EQS. 
SEMI-TRANSF  VAR,"  TBLES(  AFTER  SWITCHING). 


IN 


/BLCO/  NX,NXT,NP,NPT,NTR,IT,INVRS,NS,IP 
/BLC1/  "(101,2),U(101,2),V(101,2),W(101,2),B(101,2) 
/BLC6/  S1(101)SS2(101))S3(101),S4(101),S5(101),S6(1 
S7(101),S8(101),R1(101),R2(101),R3(101) 
COMMON  /BLC7/  C( 100 , 100) ,D( 100) ,DB( 100) ,DBP( 100) ,UE0( 100) , 


COMMON 
COMMON 
COMMON 


INT08330 
INT08340 
INT08350 
INT08360 
INT08370 
INT08380 
INT08390 
INT08400 
INT08410 
INT08420 
INT08430 
INT08440 
INT08450 
INT0S460 
INT08470 
INT08480 
INT08490 
INT08500 
INT08510 
INT08520 
INT08530 
INT08540 
)  /  DX2INT08550 
INT08560 
INT08570 
INT08580 
INT08590 
INT08600 
)  /  DX2  INT08610 
INT08620 
INT08630 
INT08640 
INT08650 
INT08660 
INT08670 
INT08680 
INT08690 
INT08700 
INT08710 
INT08720 
INT08730 
INT08740 
INT08750 
INT08760 
INT08770 
INT08780 
INT08790 
INT08800 
INT08810 
INT08820 
INT08830 
INT08840 
INT08850 
01),  INT08860 
,R4(101)INT08870 
GI      INT08880 
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COMMON  /GRD  /  ETA( 101) ,DETA( 101) ,A( 101) 

COMMON  /GTY  /  X( 101) ,UE( 100) ,P1( 100) ,P2( 100) ,CEL,CELH 


c 
c 
c 

INT089] 
INT089I 
INT089} 
INT089^ 

P1H  =  0.5  *  Pl(NX) 

DO  100  J=  2,NP 

INT089! 

FLARE   =1.0 

INT089i 

FB     =  0.5*(F(J3 

2)  +  F(J-1,2)) 

INT089 

UB     =  0.5*(U(J, 

2)  +  U(J-1,2)) 

INT089 

FVB    =  0.  5*(F(J, 

2)*V(J,2)+F(J-1,2)*V(J-1,2)) 

INT089S 

IF(UB  .LT.  0.0)  FLARE  =  0.0 

INT090( 

VB     =  0.  5*(V(J3 

2)  +  V(J-1,2)) 

INT090 

USB    =  0.5*(U(J3 

2)**2  +  U(J-1,2)**2) 

INT090: 

WSB    =  0.  5*(W(J3 

2)**2  +  W(J-1,2)**2) 

INT0903 

DERBV  =(B(J,2)*V(J,2)  -  B( J-l ,2)*V( J-1,2) )/DETA( J-l) 

INT090- 

FB4    =  0.  5*(F(J3 

1)  +  F(J-1,1)) 

INT090! 

VB4    =  0.  5*(V(J3 

l)  +  v(J-i,i)) 

INT090* 

USB4   =  0.  5*(U(J3 

1)**2  +  U(J-1,1)**2) 

INT090] 

WSB4   =  0. 5*(W(J, 

1)**2  +  W(J-1S1)**2) 

INT0901 

FVB4   =  0.  5*(F(J3 

1)*V(J,1)+F(J-1,1)*V(J-1,1)) 

INT0905 

DERBV4  =  (B(J,1)*V(J,1)  -  B( J-l , 1)*V( J-l, 1) )/DETA( J-i: 

INT091C 

S1(J)   =  CELH*(FB 

-  FB4)  +  P1H*F(J,2)  +  B(J,2)/DETA(J- 

■1) 

INT091] 

S2(J)   =  CELH*(FB 

-  FB4)  +  P1H*F(J-1,2)  -  B(J-1,2)/DETA(J-1) 

INT0912 

S3(J)   =  CELH*(VB 

+  VB4)  +  P1H*V(J,2) 

INT091J 

S4(J)   =  CELH*(VB 

+  VB4)  +  P1H*V(J-1,2) 

INT09U 

S5(J)   =  -CEL*FLARE*U(J,2) 

INT091! 

S6(J)  =  -CEL*FLARE*U(J-1,2) 

INT09K 

S7(J)   =  CEL--nv7(J,2) 

INT0911 

S8(J)   =  CEL*W(J-] 

L,2) 

INT091J 

c 

INT091! 

CRB    =  -DERBV4  +  CEL*WSB4  -  CEL*FLARE*USB4  -  Pl(NX)' 

'FVB4 

INT092C 

R2(J)   =  CRB  -  (DERBV  -  CEL*FLARE*USB  +  CEL*( VB+VB4)- 

KFB-FB4)  + 

INT092] 

+ 

CEL*WSB  +  P1(NX)*FVB) 

INT0922 

R1(J)  =  F(J-1,2) 

-  F(J,2)  +  DETA(J-1)*UB 

INT092: 

R3(J-1)=  U(J-1,2) 

-  U(J,2)  +  DETA(J-1)-'VB 

INT092* 

R4(J-1)=  W(J-1,2) 

-  W(J32) 

INT092! 

100 

CONTINUE 

INT092< 

c 

INT092; 

c 

BOUNDARY  CONDITIONS 

INT092J 

c 

Rl(l)   =  0.  0 
R2(l)   =  0.  0 
R4(NP)  =0.0 

INT092< 
INT093( 
INT093: 
INT0932 

IF(NX  . GE.  INVRS) 

GO  TO  120 

INT093: 

GAMMA  1  =  0.0 

INT093^ 

GAMMA2  =  1.  0 

INT093.' 

R3(NP)  =  0.  0 

INT093< 

RETURN 

INT093" 

120 

CONTINUE 

INT093* 

CII    =  C(NX,NX) 

*  SQRT(X(NX)) 

INT093< 

GAMMA  1  =1.0 

INT094( 

GAMMA2  =  (1.0  -  CII*ETA(NP))/CII 

INT094. 

R3(NP)  =  (GI  +  CII*(ETA(NP)*W(NP,2)  -  F(NP,2))  -W(NP,: 

0)/GII 

INT0942 

RETURN 


SO 


c 
c 
c 

c 
c 
c 
c 


c 
c 
c 

c 


END 

DATA  SET  KCBCCOEFTR  AT  LEVEL  001  AS  OF  08/24/84 
DATA  SET  KCBCCOEFTR  AT  LEVEL  001  AS  OF  08/24/84 
DATA  SET  KCBCCOEFTR  AT  LEVEL  004  AS  OF  02/21/84 

SUBROUTINE  COEFTR 

CALCULATE  COEFFS.  OF  B.  L.  FINITE-DIFFERENCE  EQS. 
IN  TRANSFORMED  VARIABLES(  BEFORE  SWITCHING). 

COMMON  /BLCO/  NX,NXT,NP,NPT,NTR, IT, INVRS ,NS , IP 
COMMON  /BLC1/  F( 101,2) ,U( 101,2) ,V( 101,2) ,W( 101,2) ,B( 101,2) 
COMMON  /BLC6/  Sl(  101) ,S2( 101) ,S3( 101) ,S4( 101) ,S5( 101) ,S6( 101) , 
+  S7(101),S8(101),R1(101),R2(101),R3(101),R4( 

COMMON  /GRD  /  ETA( 101) ,DETA( 101) , A( 101) 
COMMON  /GTY  /  X( 101) ,UE( 100) ,P1( 100) ,P2( 100) ,CEL,CELH 


DO  100  J=  2,NP 

FB     =  0.  5*(F(J,2)  +  F(J-1,2)) 

UB     =  0.5*(U(J,2)  +  U(J-1,2)) 

VB     =  0. 5*(V(J,2)  +  V(J-1,2)) 

USB    =  0. 5*(U(J,2)**2  +  U(J-1,2)**2) 

DERBV  =  (B(J,2)*V(J,2)  -  B( J-l , 2)*V( J-l , 2) )/DETA( J-l) 

FVB  =  0. 5*(F(J,2)*V(J,2)  +  F( J-1,2)*V( J-l, 2) ) 

FVB4  =  0.5*(F(J,1)*V(J,1)  +  F(J-1,1)*V(J-1,1)) 

FB4  =  0.5*(F(J,1)  +  F(J-1,1)) 

VB4  =  0.5*(V(J,1)  +  V(J-1,1)) 

USB4   =  0.  5*(U(J,1)**2  +  U(J-1,1)**2) 

DERBV4  =  (B(J,1)*V(J,1)  -  B( J-l , 1)*V( J-l , 1) )/DETA( J-l) 


S1(J)  =  CELH*(FB-FB4)  + 

S2(J)  =  CELH*(FB-FB4)  + 
+  DETA(J-l) 

S3(J)  =  CELH*(VB  +  VB4) 

S4(J)  =  CELH*(VB  +  VB4) 

S5(J)  =  -(CEL+P2(NX))*U(J,2) 

S6(J)  =  -(CEL+P2(NX))*U(J-1,2) 


0. 5*P1(NX)*F(J,2)  +  B(J,2)/DETA(J-1) 
0.5*P1(NX)*F(J-1,2)  -B(J-1,2)/ 


0. 
0. 


5*P1(NX)*V(J,2) 

5*P1(NX)*V(J-1,2) 


P2(NX-1)*USB4  +  P2(NX-1) 


CLB  =  DERBV4  +  P1(NX-1)*FVB4 

CRB  =  -CLB  -  CEL--USB4  -  P2(NX) 

R2(J)  =  CRB  -  (DERBV  +  P1(NX)*FVB-  (CEL+P2(NX))*USB 
+  (VB  +  VB4)*(FB  -  FB4)) 


+  CEL* 


R1(J)=  F(J-1,2)  -  F(J,2)  + 
R3(J-1)=  U(J-1,2)  -  U(J,2) 


DETA(J-1)*UB 
+  DETA(J-1)*VB 


100 

CONTINUE 

Rl(l)   =  0.  0 

R2(l)   =  0.  0 

R3(NP)  =0.0 

RETURN 

END 

c 

DATA  SET  KCBCCOMPBL  AT  LEVEL  001  AS  OF  08/24/84 

c 

DATA  SET  KCBCCOMPBL  AT  LEVEL  001  AS  OF  08/24/84 

INT09450 
INT09460 
INT09470 
INT09480 
INT09490 
INT09500 
INT09510 
INT09520 
INT09530 
INT09540 
INT09550 
INT09560 
101)INT09570 
INT09580 
INT09590 
INT09600 
INT09610 
INT09620 
INT09630 
INT09640 
INT09650 
INT09660 
INT09670 
INT09680 
INT09690 
INT09700 
INT09710 
INT09720 
INT09730 
INT09740 
INT09750 
INT09  760 
INT09  770 
INT09780 
INT09790 
INT09800 
INT09810 
INT09820 
INT09830 
INT09840 
INT09850 
INT09860 
INT09870 
INT09880 
INT09890 
INT09900 
INT09910 
INT09920 
INT0993O 
INT09940 
INT09950 
INT09960 
INT09970 
INT09980 
INT09990 
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c 
c 


c 

c  • 
c 

90 


110 

112 

120 
122 

130 


160 

170 

180 

182 

190 

200 

220 

210 

C 

C 

C 


DATA  SET  KCBCCOMPBL  AT  LEVEL  010  AS  OF  08/24/84 

SUBROUTINE  COMPBL( CASEID,XP, YP , VT, S ,DLSP,DLS ,DBPP,NBL, ITRI ,XCTRI , 
+  RN,NT,ISWPT) 

COMMON  /BLCO/  NX,NXT,NP ,NPT,NTR, IT, TNVRS ,NS , IP 

COMMON  /BLC7/  C( 100 , 100) ,D( 100) ,DB( 100) ,DBP( 100) ,UE0( 100) ,GI 

COMMON/ EDDY 1 /  RL , RX , SQRX , RXNTR , GMTR , GMTRS (  100) 
+  , ALFAS( 100) ,FFS( 100) ,RTS( 100) , IEDY,NXSPT 

COMMON  /GTY  /  X( 101) ,UE( 100) ,P1( 100) ,P2( 100) ,CEL,CELH 

COMMON  /BLIN/  TITLE( 20) ,XC( 100) ,YC( 100) , ISG( 100) ,DELS( 100) , 
+  XCTR,XTR,ISTRP,ICYCLE,ICYTL,XCTRS(2),TRFIND(2) 

COMMON  /BLOW/  VN(IOO) 

COMMON  /I SURF/  ISF 

C0MM0N/PL0T/NVP(2),NXVP(20,2),ICC 


DIMENSION 


DIMENSION 


+ 


XP( 100) ,DLSP( 100) ,YP( 100) ,VT( 100) ,S( 100) , 
DBPP( 100) ,DLS( 100) ,CASEID(20) 
XIN(100,2),YIN(100,2),SI(100,2),VIN(100,2), 

DIN(100,2),DELSTR(100,2),DD(100,2),DDD(100,2) 
XB(101),D1(101),D2(101),D3(101),IEND(2) 
NBL(2),ITRI(2),XCTRI(2) 


DIMENSION 

DIMENSION 

LOGICAL  TRFIND 

CHARACTER  *  4  SURF(4) ,STITLE( 2) , SURFID(4) 


DATA     SURF   /  *     ! , ' R  SU ' , ' RFAC ' , ' E 
DATA     STITLE  /  ' UPPE ' , ' LOWE '  / 


+ 


FORMAT  (  1H1 
FORMAT  (  1H0 
FORMAT  (  1H0 


FORMAT 
FORMAT 
FORMAT 


(  1H 
(  1H 
(  1H0 


+ 


140    FORMAT  (  1H0 

+ 

+ 
150   FORMAT  (  1H0 


+ 


FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 


(  1H0 

(  1H0 

(  1H0 

(  1H0 

(  1H0 

(  1H0 
(1H0, 

(  1H0 


,5X,' COMPUTING  BOUNDARY  LAYER  USING  HILBERT' , 

'  INTEGRAL',/) 
,6X,'l' ,9X,'XP' ,13X,'YP' ,13X,'S' ,14X,'VT' ,13X, 

'DBP'  /  ) 
,6X,'r ,4X,'II' ,3X,'IK' ,7X,'XIN' ,12X,'YIN' , 

13X,'Sl' ,12X,'VIN' ,12X,'DIN'  /  ) 
,5X,I3,5E15. 6  ) 
,3X,3(2X,I3),5E15. 6  ) 
,5X,' STAGNATION  POINT  IS  FOUND  BETWEEN  POINT  NO.  ', 

13,'  AND  POINT  NO.  ' ,13  /  ) 
,5X,' INTERPOLATED  STAGNATION  POINT  VALUES'  / 

1H0,5X,'S  =  ' ,E13. 6,2X,'XP  =  ' ,E13. 6 , 2X, ' YP  =  ', 

E13.  6,2X,'DBP  =  '  ,E13.  6,2X,'VT  =  0.0'  /  ) 
,5X, 'TOTAL  NUMBER  OF  UPPER  SURFACE  POINTS  =  ',15, 

2X,'AND  AT  LOWER  SURFACE  =  ',15  /  ) 
,5X, 'UPPER  SURFACE  DATA'  ) 
,5X,' LOWER  SURFACE  DATA'  ) 
,5X, 'UPPER  SURFACE  CALCULATIONS'  /  ) 
,5X, 'LOWER  SURFACE  CALCULATIONS'  /  ) 
,5X, 'RESULTS  OF  POINT  REDISTRIBUTION'  /  ) 
,5X, 'TABLE  OF  DELTA-STARS'  /  (  1H  ,5X,8E15.6  )  ) 
5X, 'TABLE  OF  BLOWING-VEL.  '  /  ( 1H  ,5X,8E15.6)  ) 
,5X,'N0  CHANGE  OF  SIGN  ON  VT.  CANNOT  FIND  STAG.  PT. ' 


INTlOOOl 
INT100H 
INT1002I 
INT1003! 
INT1004 
INT1005 
INT1006 
INT1007I 
INT1008I 
INT1009I 
INT1010 
INT1011 
INT1012 
INT1013 
INT1014 
INT1015 
INT1016 
INT1017J 
INT1018I' 
INT1019I 
INT1020( 
INT102K 
INT1022( 
INT1023( 
INT1024( 
INT1025( 
INT1026( 
INT1027( 
INT1028( 
INT1029( 
INT10301 
INT103K 
INT1032( 
INT1033( 
INT1034( 
INT1035( 
INT1036( 
INT1037( 
INT1038( 
INT1039( 
INT1040( 
INT104K 
INT1042( 
INT1043( 
INT1044( 
INT1045( 
INT1046( 
INT1047( 
INT1048( 
INT1049( 
INT1050( 
INT105K 
)INT1052( 
INT1053I 
INT1054< 
INT1055( 
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c 
c 

c 

c 


230 


300 
C 
C 
C 

c 
c 
c 

c 

C500 
C 

c 
c 


600 


READ  ONE  STRIP  INPUT  DATA  FROM  UNIT  NO.  1.   THE  ORDER  IS 
FROM  THE  LOWER  SURFACE  T.  E.  TO  THE  UPPER  SURFACE  T. E. 

DO  230  I  =  1,20 

TITLE(I)  =  CASEID(I) 

CONTINUE 

RL     =  RN 

DO  300  I  =  1,NT 

DBP(I)  =  DBPP.I) 

CONTINUE 

PRINT  THE  INPUT  DATA. 

P2(l)   =  1.  0 

WRITE  (  6,90  ) 

WRITE  (  6,110  ) 

DO  500  I  =  1,NT 

WRITE  (  6,120  )  I,XP(I),YP(I),S(I),VT(I),DBP(I) 

CONTINUE 

FIND  STAGNATION  POINT 

DO  600  I  =  2, NT 

VPROD  =  VT(I)  *  VT(I-l) 

IF  (  VPROD  . GT.  0.  )  GO  TO  600 

IS     =1 

IS1    =  IS  -  1 

GO  TO  700 

CONTINUE 

WRITE  (  6,210  ) 

STOP  1 


c 
c 
c 
c 
c 

INTERPOLATE  S  AT  VT  =  0. 

WRITE 

(  6,130  )  IS1,IS 

700 

DS 
DV 

SS 

DBB 

DX 

DY 

DSl 

DBS 

XPS 

YPS 

=  S(IS)  -  S(IS1) 

=  VT(IS)  -  VT(ISl) 

=  S(IS)  -  VT(IS)  *  (  DS/DV  ) 

=  DBP(IS)  -  DBP(ISl) 

=  XP(IS)  -  XP(ISl) 

=  YP(IS)  -  YP(ISl) 

=  SS  -  S(IS) 

=  DBP(IS)  +  DSl  *  (  DBB/DS  ) 

=  XP(IS)   +  DSl  *  (  DX  /DS  ) 

=  YP(IS)   +  DSl  *  (  DY  /DS  ) 

C 

c 
c 

WRITE 

(  6,140  )  SS, XPS, YPS, DBS 

IU  IS 

THE  TOTAL  UPPER  SURFACE  POINTS. 

+  STAG. 

PT 

c 

c 

IL  IS 

THE  TOTAL  LOWER  SURFACE  POINTS 

+  STAG. 

PT 

IU 

=  NT  -  IS  +  2 

IL 

=  IS 

c 

c 

WRITE 

(  6,150  )  IU,IL 

INT10560 
INT10570 
INT10580 
INT10590 
INT10600 
INT10610 
INT10620 
INT10630 
I NT 10 640 
INT10650 
INT10660 
INT10670 
INT10680 
INT10690 
INT10700 
INT10710 
INT10720 
INT10730 
INT10740 
INT10750 
INT10760 
INT10770 
INT10780 
INT10790 
INT10800 
INT10810 
INT10820 
INT10830 
INT10840 
INT 10 850 
INT10860 
INT10870 
INT10880 
INT10890 
INT10900 
INT10910 
INT10920 
INT10930 
INT10940 
INT10950 
INT10960 
INT10970 
INT10980 
INT10990 
INT11000 
INT11010 
INT11020 
INT11030 
INT11040 
INT11050 
INT11060 
INT11070 
INT11080 
INT11090 
INT11100 
INT11110 
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c 
c 


c 
c 
c 
c 

800 


900 


C 

1000 


1100 

c 

1200 

C 

C 

c 
c 


GROUP  THE  DATA  FOR  EACH  SURFACE.  FIRST,  UPPER. 

DO  1200  L  =  1,2 

GO  TO  (  800,900  ),  L 

L  =  1  IS  UPPER  SURFACE 
L  =  2  IS  LOWER  SURFACE 

Ml     =  IS 
M2     =  NT 
IEND(L)  =  IU 
GO  TO  1000 
Ml     =  1 
M2     =  IL-1 
IEND(L)  =  IL 

I  =  1 

XIN(I,L)  =  XPS 

YIN(I,L)  =  YPS 

SI  (I,L)  =  0. 

DIN(I,L)  ~  DBS 

VIN(I, L)  =  0. 

IF  (  IP  .GE.  1  )  THEN 

IF  (  L  .EQ.  1  )  WRITE  (  6,160  ) 

IF  (  L  . EQ.  2  )  WRITE  (  6,170  ) 

WRITE  (  6,112  ) 

WRITE  (  6,122  )  I,I,I,XIN(1,L),YIN(1,L),SI(1,L),VIN(1,L), 
+  DIN(1,L) 

END  IF 

DO  1100  II  =  M1,M2 

I      =1  +  1 

IK     =  II 

IF  (  L  .EQ.  2  )  IK  =  IL  -  II 

XIN(I,L)  =  XP(IK) 

YIN(I,L)  =  YP(IK) 

SI  (I,L)  =  S(IK)-SS 

IF  (  L  .EQ.  2  )  SI(I,L)  =  SS-S(IK) 

VIN(I, L)  =  ABS(VT(IK)) 

DIN(I.L)  =  DBP(IK) 

IF(IP  .GE.  1)V.TRITE  (  6,122  )  I ,  II ,  IK,XIN(  I  ,L)  ,YIN(  I  ,L)  ,SI(  I  ,L)  , 
+  VINCI,L),DIN(I,L) 

CONTINUE 

CONTINUE 

RE-DISTRIBUTE  POINTS  ON  EACH  SURFACE  TO  A  DENSER  NUMBER. 

WRITE  C  6,90  ) 
DO  2000  ISF  =  1,2 
NN     =  IEND(ISF) 
ITR    =  ITRI(ISF) 
NXT    =  NBLCISF) 
XCTR   =  XCTRI(ISF) 
SURF   CI)  =  STITLE(ISF) 
ICC    =  1 
DO  1220  J  =  1,4 


INT1112f 
INT1113C 
INT1114C 
INT1115C 
INT1116C 
INT1117C 
INT1118C 
INT1119C 
INT1120C 
INT1121C 
INT1122A 
INT1123C 
INT1124C 
INT1125C 
INT1126C 
INT1127C 
INT1128I 
INT1129J 
INT1130d 
INT1131C 
INT1132C 
INT1133C 
INT1134d 
INT1135d 
INT1136C 
INT1137C 
INT1138d 
INT1139d 
INT1140d 
INT1141C 
INT1142C 
INT1143C 
INT1144J 
INT1145C 
INT1146C 
INT1147C 
INT1148Q 
INT1149C 
INT1150C 
INT1151C 
INT1152C 
INT1153C 
INT1154C 
INT1155I 
INT1156C 
INT1157C 
INT1158C 
INT1159d 
INT1160C 
INT1161C 
INT1162C 
INT1163d 
INT1164J 
INT1165d 
INT11660 
INT1167C 
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SURFID(J)  =  SURF(J) 

INT11680 

1220 

CONTINUE 

INT11690 

C 

IF  (  ISF  . EQ.  1  )  WRITE  (  6,180  ) 

INT11700 

c 

IF  (  ISF  .EQ.  2  )  WRITE  (  6,182  ) 

INT11710 

SCALE   =  SI(NN,ISF) 

INT11720 

c 

INT11730 

CALL  BLGRID  (  NXT,XB,D1,D2  ) 

INT11740 

c 

INT11750 

DO  1300  I  =  1,NXT 

INT11760 

1300 

X  (I)   =  XB(I)  *  SCALE 

INT117  70 

C 

INT11780 

C 

INTERPOLATE  S,VT,X,Y,D  INTO  THE  NEW  DISTRIBUTION. 

INT11790 

C 

INT11800 

CALL  SMFIT  (  1 ,NN, SI( 1 , ISF) , VIN( 1 , ISF) ,D1 , 2  ) 

INT11810 

CALL  SMFIT  (  1 ,NN, SI( 1 , ISF) ,DIN( 1 , ISF) ,D1 ,2  ) 

INT11820 

CALL  DIFF3  (  NN,SI( 1 , ISF) , VIN( 1 , ISF) ,D1 ,D2 ,D3 , 0  ) 

INT11830 

CALL  INTRP3(  NN,SI( 1 , ISF) , VIN( 1 , ISF) ,D1 ,D2 ,D3 ,NXT,X,UE  ) 

INT11840 

CALL  DIFF3  (  NN,SI( 1 , ISF) ,DIN( 1 , ISF) ,D1 ,D2 ,D3 , 0  ) 

INT11850 

CALL  INTRP3(  NN, SI( 1 , ISF) ,DIN( 1 , ISF) ,D1 ,D2 ,D3,NXT,X,DB  ) 

INT11860 

CALL  DIFF3  (  NN, SI( 1 , ISF) ,XIN( 1 , ISF) ,D1 ,D2 ,D3 , 0  ) 

INT11870 

CALL  INTRP3(  NN,SI( 1 , ISF) ,XIN( 1 , ISF) ,D1 ,D2 ,D3 ,NXT,X,XC  ) 

INT11880 

CALL  DIFF3  (  NN,SI( 1 , ISF) , YIN( 1 , ISF) ,D1 ,D2,D3 , 0  ) 

INT11890 

CALL  INTRP3(  NN,SI( 1 , ISF) , YIN( 1 , ISF) ,D1 ,D2 ,D3 ,NXT,X,YC  ) 

INT11900 

IF  (  IP  . GE.  1  )  THEN 

INT11910 

WRITE  (  6,190  ) 

INT11920 

WRITE  (  6,110  ) 

INT11930 

DO  1320  I  =  l.NXT 

INT11940 

WRITE  (  6,120  )  I,XC(I),YC(I),X(I),UE(I),DB(I) 

INT11950 

1320 

CONTINUE 

INT11960 

END  IF 

INT11970 

C 

INT11980 

C 

INPUT  TO  THE  B. L.  PROGRAM  X,UE,DB ,XC , YC  ARE  NOW  DEFINED. 

INT11990 

C 

INT12000 

DO  1350  I  =  1,NXT 

INT12010 

DBP(I)  =  DB(I) 

INT12020 

D(I)    =  DB(I) 

INT12030 

1350 

CONTINUE 

I NT 12 040 

C 

INT12050 

CALL  BL2D(  ITR, ISWPT, SURFID) 

INT12060 

C 

INT12070 

CALL  DIFF3  (  NXT,X,DELS ,D1 ,D2 ,D3 , 0  ) 

INT12080 

CALL  INTRP3(  NXT, X, DELS ,D1 ,D2 ,D3 ,NN,SI( 1 , ISF) ,DELSTR( 1 , ISF)  ) 

INT12090 

CALL  DIFF3(NXT,X, 0,01,02,03,0) 

INT12100 

CALL  INTRP3(NXT,X,D,D1,D2,D3,NN,SI(1,ISF),DD(1,ISF)) 

INT12110 

CALL  DIFF3(NN,SI(1,ISF),DD(1,ISF),DDD(1,ISF),D2,D3,0) 

INT12120 

TRFIND(ISF)  =  .FALSE. 

INT12130 

IF(ITR  .EQ.  3  .AND.  NTR  .  LE.  NXT)  THEN 

INT12140 

XCTRS(ISF)  =  XCTR 

INT12150 

TRFIND(ISF)  =  .TRUE. 

INT12160 

END  IF 

INT12170 

C 

INT12180 

2000 

CONTINUE 

INT12190 

C 

INT12200 

C 

PUT  THE  TWO  SURFACES  DELTA- STARS  BACK  TO  ONE  STRIP 

INT12210 

C 

INT12220 

DELSTR(1,2)  =  0. 5*(DELSTR( 2 , 1)+DELSTR( 2 , 2) ) 

INT12230 
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DELSTRC 1 

DD(1,2) 

DD(1,1) 

DDD(1,2) 

DDD(1,1) 

IL  =  IEND(2) 

I      =  IL 

DO  2100  II  =  2 

I      =  1-1 

DLS(I)  =  DELSTRC 11,2) 

VN(I)  =  DDD(II,2)/SQRT(RL) 

DBPP(I)  =  DD(II,2) 


1)  =  DELSTRC 1,2) 

=0.5  *  (DD(2,1)  +  DD(2,2)) 

=  DD(1,2) 

=0.5  *  (DDD(2,1)  +  DDD(2,2))  /SQRT(RL) 

=  DDD(1,2) 


IL 


2100 

CONTINUE 

C 

I               =   IL-1 

IU            =   IEND(l) 

DO   2200   11=2, IU 

I               =1   +   1 

2200 
C 
C 
C 


C 
C 
C 

c 
c 
c 


c 
c 
c 


DLS(I)  =  DELSTR( I 1,1) 

VN(I)  =  DDD(II,1)/SQRT(RL) 

DBPP(I)  =  DD(II,1) 

CONTINUE 

WRITE  (  6,200  )  (  DLS(I),  1=1, NT  ) 

WRITE  (  6,220)  (VN(I)  ,  1=1, NT) 

RETURN 

END 

DATA  SET  KCBCCOMPGI  AT  LEVEL  001  AS  OF  08/24/84 
DATA  SET  KCBCCOMPGI  AT  LEVEL  001  AS  OF  08/24/84 
DATA  SET  KCBCCOMPGI  AT  LEVEL  003  AS  OF  08/24/84 

SUBROUTINE  COMPGI 

CALCULATE  GI 

COMMON  /BLCO/  NX,NXT,NP,NPT,NTR, IT, INVRS,NS , IP 

COMMON  /BLC7/  C( 100 , 100) ,D( 100) ,DB( 100) ,DBP( 100) ,UE0( 100) ,GI 


180 
C 


260 
C 


SUM2 

Nl 

N2 

DO  130  K 

SUM2 
CONTINUE 

Nl 

N2 

SUM1 

DO  260  K 

SUM1 

CONTINUE 

GI 

RETURN 


0.  0 

NX  -  1 

NXT 

N1,N2 

SUM2  +  C(K,NX)*  (D(K) 


2 

NX-1 

0. 

N1,N2 

SUM1  +  C(K,NX)*  (D(K) 

UEO(NX)  +  SUM1  +  SUM2 


DBP(K)) 


DBP(K)) 

C(NX,NX)  *  DBP(NX) 


INT1224C 

INT1225C 

INT1226C 

INT1227C 

INT1228J 

INT1229C 

INT1230C 

INT1231C 

INT1232C 

INT1233C 

INT1234C 

INT1235Q 

INT12360 

INT12370 

INT12380 

INT12390 

INT12400 

INT12410 

INT12420 

INT1243 

INT1244 

INT1245 

INT1246 

INT1247 

INT1248 

INT1249 

INT12500 

INT12510 

INT12520 

INT12530 

INT12540 

INT12550 

INT12560 

INT12570 

INT1258 

INT1259 

INT12600 

INT12610 

INT12620 

INT12630 

INT12640 

INT12650 

INT12660 

INT12670 

INT12680 

INT12690 

INT12700 

INT12710 

INT12720 

INT12730 

INT12740 

INT12750 

INT12760 

INT12770 

INT12780 
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86 


c 
c 
c 

c 
c 
c 

c 
c 
c 
c 
c 
c 


c 

10 


c 
c 
c 
11 


c 
c 
c 

12 


END 

DATA  SET  KCBCDIFF3  AT  LEVEL  001  AS  OF  08/24/84 
DATA  SET  KCBCDIFF3  AT  LEVEL  001  AS  OF  08/24/84 
DATA  SET  KCBCDIFF3   AT  LEVEL  002  AS  OF  04/05/84 

SUBROUTINE  DIFF3  (N,X,F,FP,FPP,FPPP, IEND) 

DETERMINES  THE  DERIVATIVE  OF  THE  INPUT  FUNCTION  AT  THE  INPUT  PTS. 
DIMENSION  X( 101) ,F( 101) ,FP( 101) ,FPP( 101) ,FPPP( 101) 


FIRST  DERIVATIVES  USING  WEIGHTED  ANGLES 

N1=N-1 
DX=X(2)-X(1) 

DF=F(2)-F(1) 
ANG2=  ATAN2(DF,DX) 
DL2=DX 

DO  10  1=2, Nl 

ANG1=ANG2 

DL1=DL2 

11=1+1 

DX=X(I1)-X(I) 

DF=F(I1)-F(I) 

ANG2=  ATAN2(DF,DX) 

DL2=DX 

ANG=(  DL2*ANG  1+DL1*ANG2  )  /  (  DL1+DL2  ) 

FP(I)=  TAN(ANG) 

IF  (I.NE.  2)  GO  TO  10 
ANGI  =  ANG 
ANG1I  =  ANGI 
DLI  =  DL1 

CONTINUE 

ANGF  =  ANG 

ANG2F  =  ANG2 

DLF  =  DL2 

IEND1  =  IEND  +  1 

GO  TO  (11,12,13),  IEND1 

IEND  =  0  ,  EXTRAPOLATE  FOR  END  VALUES 

FP(1)  =  2.*(F(2)-F(1  ))/DLI  -  FP(2) 
FP(N)  =  2.*(F:N)-F(N1))/DLF  -  FP(N1) 
GO  TO  20 

IEND  =  1,  DERIVATIVES  ARE  CONTINUOUS  ACROSS  ENDS 

ANG  -  (DLI*ANG2F  +  DLF*ANG1I)  /  (DLI  +  DLF) 
FP(1)  =  TAN(ANG) 
FP(N)  =  FP(1) 


INT12790 
INT12800 
INT12810 
INT12820 
INT12830 
INT12840 
INT12850 
INT12860 
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INT12880 
INT12890 
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INT12950 
INT12960 
INT12970 
INT12980 
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INT13010 
INT13020 
INT13030 
INT13040 
INT13050 
INT13060 
INT13070 
INT13080 
INT13090 
INT13100 
INT13110 
INT13120 
INT13130 
INT13140 
INT13150 
INT13160 
INT13170 
INT13180 
INT13190 
INT13200 
INT13210 
INT13220 
INT13230 
INT13240 
INT13250 
INT13260 
INT13270 
INT13280 
INT13290 
INT13300 
INT13310 
INT13320 
INT13330 
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c 
c 
c 

13 


c 
c 
c 

20 


30 


C 

c 
c 

c 
c 
c 


c 
c 
c 


GO  TO  20 

IEND  =  2,  IF  FIRST  DERIVATIVE  . LT.  0.0  RESET  TO  ZERO 

CONTINUE 

FP(1)  =  2.*(F(2)-F(1  ))/DLI  -  FP(2) 
IF  (FP  (1)  .  LT.  0.0)  FP  (1)  =  0.0 
FP(N')  =  2.*(F(N)-F(N1))/DLF  -  FP(N1) 

SECOND  +  THIRD  DERIVATIVES  USING  CUBIC  FITS 


10 


30 


DO 

II 

12 

DX1 

DX2 

DF1 

DF2 

FPPP(I)= 

FPP  (I)= 

CONTINUE 

FPPP(1)= 

FPPP(N)= 

FPP  (1)= 

FPP  (N)= 


1=2, Nl 
=  1-1 
=  1  +  1 
=  X  (II) 
=  X  (12) 
=  2.0  * 
=  2.  0  * 
3.0  * 
DF1  -  DX1 


-  X  (I) 

-  X  (I) 
((F  (II) 
((F  (12) 
(DF1 


-  F  (I))  /  DX1  - 

-  F  (I))  /  DX2  - 
DF2)  /  (DX1  -  DX2) 

FPPP  (I)  /  3.  0 


FP 
FP 


CD) 

(D) 


DXl 
DX2 


FPPP  (2) 
FPPP  (Nl) 
FPP  (2  )  +  (X 
FPP  (Nl)  +  (X 


(1) 
(N) 


X 
X 


(2  )) 

(Nl)) 


FPPP 
FPPP 


RETURN 

END 

DATA  SET  KCBCEDDY 
DATA  SET  KCBCEDDY 
DATA  SET  KCBCEDDY 

SUBROUTINE  EDDY 


AT  LEVEL 
AT  LEVEL 
AT  LEVEL 


001 
001 
003 


AS 
AS 
AS 


OF 
OF 
OF 


(2  ) 
(Nl) 


08/24/84 
08/24/84 
04/05/84 


CALCULATE  EDDY  VISCOSITY  USING  C  . S.  TWO-LAYERS  EDDY 

VISCOSITY  FORMULA 

COMMON  /BLCO/  NX,NXT,NP,NPT,NTR, IT, INVRS ,NS , IP 
COMMON  /BLC1/  F( 101 , 2) ,U( 101 , 2) , V( 101 , 2) ,W( 101 ,2) ,B( 101, 2) 
COMMON/BLC7/  C( 100 , 100) ,D( 100) ,DB( 100) ,DBP( 100) ,UE0( 100) ,GI 
COMMON/EDDY1/  RL,RX,SQRX,RXNTR,GMTR,GMTRS( 100) , 
+  ALFAS( 100) ,FFS( 100) ,RTS( 100) , IEDY,NXSPT 

COMMON  /GRD  /  ETA( 101) ,DETA( 101) ,A( 101) 
COMMON  /GTY  /  X( 101) ,UE( 100) ,P1( 100) ,P2( 100) ,CEL,CELH 
DIMENSION  FINT(lOl) 


JO=l 

UDEL=0.995*U(NP,2) 
DO  10  J=1,NP 
IF(U(J,2).LT.  UDEL)  JJ=J 
IF(U(J,2).LT.  0.0)   JO=J 
EDEL=ETA(JJ)+(ETA(JJ+1)- 
+     *(UDEL-U(JJ,2)) 
DO  15  J=1,NP 
ETADEL=ETA(J)/EDEL 


ETA(JJ))/(U(JJ+1,2)-U(JJ,2)) 


INT1334 

INT1335 

INT1336i 

INT1337' 

INT1338 

INT1339" 

INT1340I 

INT134H 

INT1342 

INT134: 

INT134^ 

I  NT  13  451 

INT134* 

INT13471 

INT1348C 

INT1349C 

INT1350C 

INT1351C 

INT1352! 

INT1353J 

INT1354C 

INT1355C 

INT1356C 

INT1357( 

INT1358J 

INT1359( 

INT13600 

INT13610 

INT13620 

INT13630 

INT13640 

INT13650 

INT13660 

INT1367( 

INT1368C 

INT1369( 

INT13700 

INT13710 

INT13720 

INT13730 

INT13740 

INT13750 

INT13760 

INT13770 

INT13780 

INT13790 

INT13800 

INT13810 

INT13820 

INT13830 

INT13840 

INT13850 

INT13860 

INT13870 

INT13880 


88 


15 


20 


30 


C 
C 

r 

35 


IF(ETADEL.  GT.  1.  0)  ETADEL=1.  0 
FINTC  J)  =  l.  0/(1.  0+5.  5*ETADEL**6) 

CALL  AME AN ( 1 , J J , ETA , F I NT , 2 ) 
RU  =  RL 

IF  (IT  . GT.  1)  GO  TO  20 
ALFAS(NX)  =  ALFAS(NX-l) 
FFS(NX)  =  FFS(NX-l) 
RTS(NX)  =  RTS(NX-l) 

GMTR   =  GMTRS(NX) 
IF  (NX  .LE.  NS)  RU  =  RL  *  UE(NX) 
RL2    =  SQRT(RU  *  X(NX)) 
RL4    =  SQRT(RL2) 
RL216  =  0.  16  *  RL2 
VMAX   =0.5  *  (V(l,2)  +  V(l,l)) 
DO  30  J=2,NP 

VM     =0.5  *  (V(J,2)+V(J,1)) 
IF(VM  .GT.  VMAX)  VMAX  =  VM 
CONTINUE 

IF  (IEDY  .EQ.  0)  GO  TO  35 
IF  (IT  .LE.  1  .OR.  GMTR  .  LT.  0.85 
I-     GO  TO  35 


OR.  NX  . LE.  NTR+3  ) 


MODIFY  ALFA  USING  SIMPSON'S  ARGUMENTS 
CALL  SMPSON 
ALFA  =  ALFAS(NX) 

EDVO   =  ALFA  *  RL2  *  GMTR  *  (U(NP,2)*ETA(NP)  -  F(NP,2)) 
DO  40  J=2,NP 
JJ     =  J 

YBA    =  RL4*SQRT(VMAX)/26.0*ETA(J) 
EL     =  1.  0 

IF(YBA  .  LT.  10.0)  EL  =  1.  0  -  EXP(  -YBA) 
EDVI   =  RL216  *  GMTR  *  (EL*ETA( J))**2  *   ABS(V(J,2)) 
IF(EDVI  .GT.  EDVO)  GO  TO  70 
B(J,2)  =  1.0  +  EDVI-FINT(J) 
IF(B(J,2)  .LT.  B(J-1,2))  B(J,2)  =  B(J-1,2) 


C 

B(J,2)  =  1.  0  +  EDVI 

40 

CONTINUE 

JM     =  2 

BM     =  B(2,2) 

DO  50  J=2,NP 

1F(BM. GT. B(J,2))  GOTO  50 

JM     =  J 

BM     =  B(J,2) 

50 

CONTINUE 

GOTO  90 

70 

DO  80  J=JJ,NP 

80 

B(J,2)  =  1.0  +  EDVO*FINT(J) 

C 

80  B(J,2)  =  1. 0  +EDVO 

c 

90 


CONTINUE 
B(l,2)  =  1.  0 

JJ     =  1 
DO  100  J=2,NP 


INT13890 
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INT13970 
INT13980 
INT13990 
INT14000 
INT14010 
INT14020 
INT14030 
INT14040 
INT14050 
INT14060 
INT14070 
INT14080 
INT14090 
INT14100 
INT14110 
INT14120 
INT14130 
INT14140 
INT14150 
INT14160 
INT14170 
INT14180 
INT14190 
INT14200 
INT14210 
INT14220 
INT14230 
INT14240 
INT14250 
INT14260 
INT14270 
INT14280 
INT14290 
INT14300 
INT14310 
INT14320 
INT14330 
INT14340 
INT14350 
INT14360 
INT14370 
INT14380 
INT14390 
INT14400 
INT14410 
INT14420 
INT14430 
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IF(U(J,2)  .LT.  0.  0)  JJ  =  J 
100    CONTINUE 

IF(JJ.  EQ.  1)  GO  TO  110 
C 

C    IN  THE  SEPARATED  REGION,  EDDY  VISCOSITY  IS  SET  EQUAL  TO 
C   THAT  IN  THE  PREVOUS  STATION  TO  AVOID  NUMERICAL  TROUBLE 

JJP3   =  JJ  +  3 

JJP3   =  MIN(JJP3,  NP) 

CALL  AMEAN(1,JJP3,ETA,B(1,2),2) 

CALL  AMEAN(1,NP,ETA,B(1,2),1) 

RETURN 

END 

DATA  SET  KCBCEXTRAP  AT  LEVEL  001  AS  OF  08/24/84 
DATA  SET  KCBCEXTRAP  AT  LEVEL  001  AS  OF  08/24/84 
DATA  SET  KCBCEXTRAP  AT  LEVEL  008  AS  OF  02/13/84 

SUBROUTINE  EXTRAP( NX , NXTE , X , Y) 

EXTRAPOLATE  DATA  USING  PARABOLIC  FORMULA 
DIMENSION  X(101),Y(101) 


110 

c 


c 
c 
c 

c 
c 


10 


c 
c 

c 


c 
c 
c 


c 
c 


=  Y(NX-2) 

=  Y(NX-l) 

=  X(NX-2) 

=  X(NX-l) 

=  X(NXTE) 

=  XI  -X3 

=  X2  -X3 

=  (Y1-Y2)/(X1**2  -  X2**2) 

=  Yl  -  BB  *  Xl**2 

I=NX,NXTE 

=  X(I)  -X3 

=  AA  +  BB  *  XI  — 2 


Yl 

Y2 

XI 

X2 

X3 

XI 

X2 

BB 

AA 

DO  10 

XI 

Y(I) 

CONTINUE 

RETURN 

END 

DATA  SET  KCBCFILLUP 

DATA  SET  KCBCFILLUP 

DATA  SET  KCBCFILLUP 
SUBROUTINE  FILLUP( INDEX) 

COMMON  /BLCO/  NX,NXT,NP,NPT,NTR, IT, INVRS ,NS , IP 
COMMON  /BLC1/  F( 101 ,2) ,U( 101 , 2) ,V( 101 ,2) ,W( 101 ,2) ,B( 101,2) 
COMMON  /GRD  /  ETA( 101) ,DETA( 101) ,A( 101) 


IF(NP.  GE.NPT  )  RETURN 
IF( INDEX. EQ. 2)  GOTO  10 

DEFINE  PROFILES  FOR  B.  L.  GROWTH 


AT 

LEVEL 

001 

AS 

OF 

08/24/84 

AT 

LEVEL 

001 

AS 

OF 

08/24/84 

AT 

LEVEL 

007 

AS 

OF 

04/05/84 

NP1 
NP 
NP 
NM 
GOTO  20 


=  NP  +  1 

=  NP  +  2 

=  MIN(NP,  NPT) 

=  NP 


INT14440 

INT14450 

INT1446G; 

INT1447CI 

INT144801 

INT14490; 

INT145O0I 

INT14510 

INT14520I 

INT14530 

INT14540 

INT14550 

INT14560 

INT14570 

INT14580 

INT14590 

INT14600 

INT14610 

INT14620 

INT14630 

INT14640 

INT14650 

INT14660 

INT14670 

INT14680 

INT14690 

INT14700 

INT14710 

INT14720 

INT14730 

INT14740 

INT14750 

INT14760 

INT14770 

INT14780 

INT14790 

INT14800 

INT14810 

INT14820 

INT14830 

INT14840 

INT14850 

INT14860 

INT1487I 

INT14880 

INT14890 

INT149O0 

INT14910 

INT14920 

INT14930 

INT14940 

INT14950 

INT14960 

INT14970 

INT14980 

INT14990 


90 


c 

FILL  UP  PROFILES  BEFORE  MOVING  TO  THE  NEXT  STATION 

INT15000 

c 

INT15010 

10 

NP1    =  NP  +  1 

INT15020 

NM     =  NPT 

INT15030 

20 

DO  30  J=NP1,NM 

I NT 15 040 

F(J,2)  =  F(J-1,2)  +  DETA(J-1)*U(J-1,2) 

INT15050 

U(J,2)  =  U(J-1,2) 

INT15060 

V(J,2)  =0.0 

INT15070 

B(J,2)  =  B(J-1,2) 

INT15080 

W(J,2)  =  W(J-i,2) 

INT15090 

30 

CONTINUE 

INT15100 

C 

INT15110 

RETURN 

INT15120 

END 

INT15130 

C 

DATA  SET  KCBCHEADER  AT  LEVEL  001  AS  OF  08/24/84 

INT15140 

C 

DATA  SET  KCBCHEADER  AT  LEVEL  001  AS  OF  08/24/84 

INT15150 

C 

DATA  SET  KCBCHEADER  AT  LEVEL  001  AS  OF  04/05/84 

INT15160 

SUBROUTINE  HEADER  (  TITLE ,SURFID , ISTRP  ) 

INT15170 

COMMON  /I SURF/  ISF 

INT15180- 

C 

INT15190 

DIMENSION     TITLE(20),   SURFID(4) 

INT15200 

C 

INT15210 

10 

FORMAT  (  1H1,20X,20A4  ) 

INT15220 

20 

FORMAT  (  1H0,15X, 'BOUNDARY  LAYER  CALCULATION  FOR  * , 

INT15230 

+         'UPPER  SURFACE  ' , 10X, ' ICYCLE=' , 15  /  16X,71(1H-)  /  ) 

INT15240 

30 

FORMAT  (  1H0,15X, 'BOUNDARY  LAYER  CALCULATION  FOR  ', 

INT15250 

+         'LOWER  SURFACE  ' , 10X, ' ICYCLE=' , 15  /  16X,71(1H-)  /  ) 

INT15260 

C 

INT15270 

p 

INT15280 
INT15290 

c 

WRITE  (  6,10  )  TITLE 

INT15300 

IF(ISF  .EQ.  1)  WRITE  (  6,20  )  ISTRP 

INT15310 

IF(ISF  .EQ.  2)  WRITE  (  6,30  )  ISTRP 

INT15320 

c 

INT15330 

RETURN 

INT15340 

END 

INT15350 

c 

DATA  SET  KCBCINPUT  AT  LEVEL  001  AS  OF  08/24/84 

INT15360 

c 

DATA  SET  KCBCINPUT  AT  LEVEL  001  AS  OF  08/24/84 

INT15370 

G 

DATA  SET  KCBCINPUT  AT  LEVEL  009  AS  OF  08/24/84 

INT15380 

SUBROUTINE  INPUT( ITR, ISWPT,SURFID) 

INT15390 

COMMON  /BLCO/  NX,NXT,NP,NPT,NTR, IT, INVRS ,NS , IP 

INT15400 

COMMON  /BLC1/  F( 101 , 2) ,U( 101 , 2) , V( 101 ,2) ,W( 101 , 2) , B( 101 , 2) 

INT15410 

COMMON  /BLC2/  DELF( 101) ,DELU( 101) ,DELV( 101) ,DELW( 101) 

INT15420 

COMMON  /BLC7/  C( 100, 100) ,D( 100) ,DB( 100) ,DBP( 100) ,UE0( 100) ,GI 

INT15430 

COMMON  /BONV/  ITMAX,EPSL,EPST,CONV 

INT15440 

COMMON/ EDDY 1/  RL ,RX , SQRX , RXNTR , GMTR , GMTRS( 100 ) , 

INT15450 

+                  ALFAS( 100) ,FFS( 100) ,RTS( 100) , IEDY,NXSPT 

INT15460 

COMMON  /GTY  /  X( 101) ,UE( 100) ,P1( 100) ,P2( 100) ,CEL,CELH 

INT15470 

COMMON  /BLIN/  TITLE( 20) ,XC( 100) , YC( 100) , ISG( 100) ,DELS( 100) , 

INT15480 

+             XCTR,XTR, ISTRP, ICYCLE , ICYTL,XCTRS( 2) ,TRFIND( 2) 

INT15490 

COMMON/TRN/  PGAMTR , OMEGA , RTHETB , RTRANB 

INT15500 

COMMON  /I SURF/  ISF 

INT15510 

DIMENSION  Dl( 100) ,D2( 100) ,D3( 100) 

INT15520 

DIMENSION  SURFID(4),XCS(100) 

INT15530 

LOGICAL  TRFIND 

INT15540 

C 

INT15550 

91 


c 
c 

UMAX 

=  15 

EPSL 

=  0.  0001 

EPST 

=  0.  01 

NPT 

=  101 

ETAE 

=  8.0 

VGP 

=  1.  14 

DETA1 

=  0.  01 

NSS 

=  NXT  /  4 

c 

SEARCH  FOR  PRESSURE  PEAK  AS  SWITCH  POINT 

UMAX  = 

UE(1) 

DO  50  ] 

[  =  2  ,  NXT 

IF  (UE(I)  . LE.  UMAX)  GO  TO  55 

UMAX  = 

UE(I) 

50 

CONTINUE 

GO  TO  60 

55 

NS  =  I 

-  4 

60 

IF  (NS 

.GT.  NSS)  NS  =  NSS 

IF  (NS 

.  LT.    3)  NS  =  3 

c 
c 
c 
c 


65 


70 


80 

C 
C 


C 
C 
95 

100 


CALCULATE  THE  PRESSURE  PARAMETERS  PI  +  P2  FOR  B. 
USING  TRANSFORMED  COORDINATES 

CALL  DIFF3  (NXT,  X,  UE ,  Dl,  D2 ,  D3 ,  0  ) 

DO  65  I  =  2, NXT 

P2(I)    =  X(I)  *  Dl( I)  /UE(I) 

P1(I)    =0.5  *  (1.0  +  P2(I)) 

CONTINUE 

Pl(l)    =  0.  5  *  (1.0  +  P2(l)) 

XCMIN   =  XC(1) 

MIN     =  1 

DO  70  1=1, NXT 

IF( XCMIN. LT. XC(I))  GOTO  70 

XCMIN   =  XC(I) 

MIN     =  I 

CONTINUE 

DO  80  I  =  1  ,  NXT 

IF  (I  . GE.  MIN)  THEN 

XCS(I)  =  XC(I) 

ISG(I)  =  1 
ELSE 

XCS(I)  =  -XCS(I) 

ISG(I)  =  -1 
END  IF 
CONTINUE 
INVRS   =  NS  +  1 

SEARCH  FOR  TRANSITION  LOCATION 
ITRP1   =  ITR  +  1 
GOTO  (150,  95,  120,  150),  ITRP1 

TRANSITON  LOCATION  IS  INPUT  (  =  XCTR) 
DO  100  1=1, NXT 
IF(XCTR. LT. XCS(I))  GOTO  105 
CONTINUE 


L.  CALCULATION 


90 
,00 
,10 


INT15560 

INT15570 

INT15580 

INT15590 

INT15600 

INT15610 

INT15620 

INT15630 

INT 15 640 

INT15650 

INT15660 

INT15670 

INT15680 

INT15690 

INT1570 

INT15710 

INT15720 

INT15730 

INT15740 

INT15750 

INT15760 

INT15770 

INT15780 

INT15790 

INT1580 

INT158 

INT15820 

INT15830 

INT15840 

INT15850 

INT15860 

INT15870 

INT15880 

INT15890 

INT15900 

INT15910 

INT1592J 

INT1593 

INT1594 

INT15950 

INT15960 

INT15970 

INT15980 

INT15990 

INT16000 

INT16010 

INT16020 

INT16030 

INT16040 

INT16050 

INT16060 

INT16070 

INT16080 

INT16090 

INT16100 

INT16110 


I 


92 


105 


C 

c 
c 

120 


NTR     =  NXT  +  1 

GOTO  200 

NTR  =1-1 

IF  (NTR.  LT.  3)  THEN 

NTR  =  3 

XCTR  =  XC(NTR) 
END  IF 
GOTO  200 

TRANSITION  LOCATION  IS  SET  AT  THE  PRESSURE  PEAK 


75 


TRANSTION  LOCATION  WILL  BE  CALCULATED  BASED  ON  MICHEL  CRITERION 


UM      =  UE(1) 
IM      =  1 

DO  75  I  =  1,NXT 

IF(UM. GT. UE(I))  GOTO  75 

IM      =  I 

UM      =  UE(IM) 

CONTINUE 

IF(IM.  LT.  4)  IM  =  4 

NTR  =  IM 

XCTR  =  XC(NTR) 

GOTO  200 
C 
C 

c 

150   NTR     =  NXT  +  1 

200   DO  210  I=1,NXT 

210   GMTRS(I)=  0.  0 

C 

C  TRANSITION  LOCATION  PROVISIONALLY  FROM  PREVIOUS  CYCLE 

C 

IF  (TRFIND(ISI)  )  THEN 

DO  211  I  =  1  ,  NXT 

XCS(I)  =  xn/T) 

IF  (I  .Li.  MIN)  XCS(I)  =  -XCS(I) 

CONTINUE 

DO  215  1=1, NXT 

IF(XCTRS(ISF)  .LE.XCS(I))  GOTO  217 

CONTINUE 

NTR  =  1-1 

XCTR  =  XCTRS(ISF) 

IF  (NTR  .LT.  3)  THEN 
NTR  =  3 
XCTR  =  XC(NTR) 

END  IF 
END  IF 


211 


215 
217 


C 
C 
C 


CALCULATE  GAMTR  DISTRIBUTION 

IF  (NTR. GT. NXT-1)  GOTO  250 

FAC  =  (XCTR-XC(NTR))/(XC(NTR+1)-XC(NTR)) 

XTR  =  X(NTR)  +  FAC*(X(NTR+1)-X(NTR)) 

UETR  =  UE(NTR)  +  FAC*(UE(NTR+1) -UE(NTR) ) 

RXNTR  =  XTR  *  UETR  *  RL 

GGFT   =  1. 0/PGAMTR*RL**2/RXNTR**l. 34 

DO  220  I=NTR+1,NXT 


INT16120 
INT16130 
INT16140 
INT16150 
INT16160 
INT16170 
INT16180 
INT16190 
INT16200 
INT16210 
INT16220 
INT16230 
INT16240 
INT16250 
INT16260 
INT16270 
INT16280 
INT16290 
INT16300 
INT16310 
INT16320 
INT16330 
INT16340 
INT16350 
INT16360 
INT16370 
INT16380 
INT16390 
INT16400 
INT16410 
INT16420 
INT16430 
INT16440 
INT16450 
INT16460 
INT16470 
INT16480 
INT16490 
INT16500 
INT16510 
INT16520 
INT16530 
INT16540 
INT16550 
INT16560 
INT16570 
INT16580 
INT16590 
INT16600 
INT16610 
INT16620 
INT16630 
INT16640 
INT16650 
INT16660 
INT16670 
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220 


230 

250 

C 

C 

C 


260 

C 

C 


ALFAS(I)  =  0. 0168 

GMTRS(I)=  1.  0 

ALFAS(NTR)  =  0. 0168 

UEINTG  =  0.0 

Ul  =  0. 5/UETR 

XI  =  XTR 

DO  230  I=NTR+1,NXT 

U2  =  0. 5/UE(I) 

X2  =  X(I) 

UEINTG  =  UEINTG+(U1+U2)*(X2-X1) 

Ul  =  U2 

XI  =  X2 

GG  =  GGFT*UEINTG*(X(I)-XTR)*UE(I)**3 

IF(GG  .GT.  10.  0)  GO  TO  250 

GMTRS(I)  =  1.  O-EXP(-GG) 

CONTINUE 

CONTINUE 

GENERATE  B.  L.  ETA  GRIDS  +  INTIAL  VELOCITY  PROFILES 

CALL  INTL(ETAE,DETA1,VGP) 
DO  260  1=1, NXT 
UEO(I)   =  UE(I) 
CONTINUE 
PRINT  OUT  INPUT  DATA 

IF  (ICYCLE  .GE.  ICYTL-1  .OR.  IP  .  GE.  0)  THEN 
CALL  HEADER(  TITLE ,SURFID, ISTRP  ) 
WRITE( 6 , 1002)  NXT, ITR, IP,NS ,NTR, ISWPT 
WRITE(6,1003)  VGP,DETA1,RL,XCTR, OMEGA, PGAMTR 

ELSE 

IF  (ISF.EQ. 1)  WRITE(6,1008)  ICYCLE 
IF  (ISF.EQ. 2)  WRITE(6,1009)  ICYCLE 

END  IF 

IF  (NTR. LT.  NXT)  THEN 

IF  (ITR. EQ. 1)  WRITE  (6,1005)  XCTR,XTR,NTR 
IF  (ITR. EQ. 2)  WRITE  (6,1006)  XCTR,XTR,NTR 
IF  (TRFIND(ISF))  WRITE( 6, 1007)  XCTR,XTR,NTR 

END  IF 

RETURN 


c 

2 

F0RMAT(20A4) 

3 

F0RMAT(10I^ 

4 

FORMAT(6F10. 0) 

1001 

FORMAT(1H1,20X, 

20A4) 

1002 

F0RMAT(1H0,10H 

NXT 

+            IH 

,10H 

+            IH 

,10H 

1003 

F0RMAT(1H0,10H 

VGP 

+            IH 

,10H 

+            IH 

,10H 

=  ,I5,7X,10H    ITR  =  ,I5,7X/ 
IP   =  ,I5,7X,10H   NS   =  ,15, 7X/ 
NTR  =  ,I5,7X,10H   ISWPT=  ,15) 
=  ,E12.4,10H   DETA1=  ,E12.  4/ 
RL   =  ,E12.4,10H   XCTR  =  ,E12.4/ 
OMEGA  =  ,E12.4,10H  PGAMTR=  ,E12. 4) 

1004  FORMAT(lH0,3X,2H  I ,6X,2HXC, 11X,2HYC, 11X,2H  X, 11X,2HUE, 11X,2HP1, 
+  11X,2HP2,/(1H  ,3X,I3,6E13.  5)) 

1005  F0RMAT(/3X,f BEGIN  OF  TRANSITION  IS  BEING  INPUT  AT  X/C  =',F8.4,4X, 


INT16680 
INT16690 
INT16700 
INT16710 
INT16720 
INT16730 
INT 16 740 
INT16750 
INT16760 
INT16770 
INT16780 
INT16790 
INT16800 
INT16810 
INT16820 
INT16830 
INT16840 
INT16850 
INT16860 
INT16870 
INT16880 
INT16890 
INT16900 
INT16910 
INT16920 
INT16930 
INT16940 
INT16950 
INT16960 
INT16970 
INT16980 
INT16990 
INT17000 
INT17010 
INT17020 
INT17030 
INT 17 040 
INT17050 
INT17060 
INT17070 
INT17080 
INT17090 
INT17100 
INT17110 
INT17120 
INT17130 
INT17140 
INT17150 
INT17160 
INT17170 
INT17180 
INT17190 
INT17200 
INT17210 
INT17220 
INT17230 
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+       'S/C  =' ,F8.4,4X,'NTR  =' ,13/) 

1006  F0RMAT(/3X, 'BEGIN  OF  TRANSITION  IS  SET  AT  PRESSURE  PEAK,  X/C  = 
+       F8.4,4X,'S/C  =' ,F8.4,4X,'NTR  =',I3/) 

1007  FORMAT(/3X, 'BEGIN  OF  TRANSITION  IS  PROVISIONALLY  TAKEN  FROM  ', 
+   'PREVIOUS  CYCLE:  X/C=' ,F8. 4 ,4X, ' S/C  =' ,F8. 4 ,4X, ' NTR  =',I3/) 

1008  FOR!'!AT(//3X,' UPPER  SURFACE  CALCULATIONS  OF  CYCLE', 13) 

1009  FORMAT(//3X,' LOWER  SURFACE  CALCULATIONS  OF  CYCLE', 13) 
END 

C  DATA  SET  KCBCINTL   AT  LEVEL  001  AS  OF  08/24/84 

C  DATA  SET  KCBCINTL   AT  LEVEL  001  AS  OF  08/24/84 

C  DATA  SET  KCBCINTL   AT  LEVEL  009  AS  OF  02/22/84 

SUBROUTINE  INTL(ETAE ,DETA1 , VGP) 
COMMON  /BLCO/  NX,NXT,NP,NPT,NTR, IT, INVRS ,NS , IP 
COMMON  /BLC1/  F( 101 , 2) ,U( 101 , 2) , V( 101 ,2) ,W( 101 , 2) ,B( 101 , 2) 
COMMON  /BLC2/  DELF( 101) ,DELU( 101) ,DELV( 101) ,DELW( 101) 
COMMON  /BLC6/  Sl( 101) ,S2( 101) ,S3( 101) ,S4(  101) ,S5( 101)  ,S6( 101) , 
+  S7(101),S8(101),R1(101),R2(101),R3(101),R4( 

COMMON  /BONV/  ITMAX,EPSL,EPST,CONV 
COMMON  /GRD  /  ETA( 101) ,DETA( 101) , A( 101) 
COMMON  /GTY  /  X( 101) ,UE( 100) ,P1( 100) ,P2( 100) ,CEL,CELH 
C 

C    

c 
c 
c 


10 
20 


30 


40 
C 

c 

80 


GENERATE  THE  GRID 

DETA(l)  =  DETA1 

IF(VGP.  LT.  1.  0)  VGP  =  1.  0 

IF((VGP-1.0)  .LE.  0.001)   GOTO  10 

NP     =  ALOG((ETAE/DETA(l))*(VGP-l. 

GO  TO  20 

NP     =  ETAE/DETA(1)  +  1.001 

IFCNP  . LE.  NPT)   GO  TO  30 

WRITE(6,  150  ) 

STOP 

ETA(l)  =  0.  0 

DO  40  J=2,NPT 

ETA(J)  =  ETA(J-l)  +  DETA(J-l) 

DETA(J)=  VGP*DETA(J-1) 

A( J)    =  0. 5*DETA(J-1) 

CONTINUE 

GENERATE  INITIAL  VELOCITY  PROFILE 
DO  90  J=1,NP 
ETAB  =  ETA(J)/ETA(NP) 
ETAB2  =  ETAB**2 


90 


120 


0)+l.  0)/ALOG(VGP)  +  1.001 


F(J,2)  = 

U(J,2)  = 

V(J,2)  = 

B(J,2)  = 

W(J,2)  = 

CONTINUE 

NX 

IT 

IT 


0.  25*ETA(NP)*ETAB2*(3.0  - 
0.5*ETAB*(3.  0  -  ETAB 2) 

-  ETAB2)/ETA(NP) 


0.5*ETAB2) 


1. 
1. 
1. 

1 
0 
IT 


'(1.0 


+  1 


IF(IT  .LT.  ITMAX)  GO  TO  130 
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INT17380 
INT17390 

101)INT17400 
INT17410 
INT17420 
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130 
C 


140 


C 
150 

C 
C 
C 

C 
C 
C 
C 
C 
C 

c 

c 

c 

c 

c 

c  - 

c 


STOP 
CONTINUE 

DO  140  J= 

FB 

UB 

VB 

USB 

DERBV  = 

FVB 

S1(J)   = 

S2(J)   = 

S3(J) 

S4(J) 

S5(J) 

S6(J) 

CRB 

R2(J) 


2,NP 
0.5*(F(J,2) 
0.5*(U(J,2) 
0.5*(V(J,2) 
0.5*(U(J,2) 


+  F(J-1,2)) 

+  U(J-1,2)) 

+  V(J-1,2)) 

«'-2  +  U(J-1,2)< 


2) 

(B(J,2)*V(J,2)  -  B(J-1,2)*V(J-1,2))/DETA(J-1) 
0.5*(F(J,2)*V(J,2)  +  F(J-1,2)*V(J-1,2)) 


0.5*P1(NX)*F(J,2)  + 
0.5*P1(NX)*F(J-1,2) 

0. 5*P1(NX)*V(J,2) 

0.5*P1(NX)*V(J-1,2) 

-P2(NX)*U(J,2) 

-P2(NX)*U(J-1,2) 

-P2(NX) 

CRB  -  (DERBV  +  P1(NX)*FVB 


B(J,2)/DETA(J-1) 

-  B(J-1,2)/DETA(J-1) 


P2(NX)*USB) 


F(J-1,2) 
U(J-1,2) 

0.  0 
0.  0 
0.  0 


F(J,2)  +  DETA(J-1)*UB 
U(J,2)  +  DETA(J-1)*VB 


GT.  EPSL   )  GO  TO  120 


R1(J)  = 
R3(J-1)= 
CONTINUE 

RKD   = 
R2(l)  = 
R3(NP)  = 
CALL  S0LV3 
IF(ABS(DELV(1)) 
CALL  FILLUP(2) 
CALL  OUTPUT(l) 

RETURN 


FORMAT(1HO,37HNP  EXCEEDED  NPT  --  PROGRAM  TERMINATED) 

END 

DATA  SET  KCBCINTRP3  AT  LEVEL  001  AS  OF  08/24/84 
DATA  SET  KCBCINTRP3  AT  LEVEL  001  AS  OF  08/24/84 
DATA  SET  KCBCINTRP3  AT  LEVEL  003  AS  OF  04/05/84 

SUBROUTINE  INTRP3  (Nl ,X1 ,F1 ,FP1 ,FPP1 ,FPPP1,N2,X2 ,F2) 

CUBIC  INTERPOLATION 

GIVEN  THE  VALUES  OF  A  FUNCTION  (Fl)  AND  ITS  DERIVATIVES 
AT  Nl  VALUES  OF  THE  INDEPENDENT  VARIABLE  (XI) 

FIND  THE  VALUES  OF  THE  FUNCTION  (F2) 

AT  N2  VALUES  OF  THE  INDEPENDENT  VARIABLE  (X2) 

X2  CAN  BE  IN  ARBITRARY  ORDER 


DIMENSION  X1(101),F1(101))FP1(101),FPP1(101),FPPP1(101),X2(101) 
+  F2(101) 

DATA  EPS  /IE -04/ 


INT17790 
INT17800 
INT17810 
INT17820 
INT17830 
INT17840 
INT17850 
INT17860 
INT17870 
INT17880 
INT17890 
INT17900 
INT17910 
INT17920 
INT17930 
INT17940 
INT17950 
INT17960 
INT17970 
INT17980 
INT17990 
INT18000 
INT18010 
INT18020 
INT18030 
INT18040 
INT18050 
INT18060 
INT18070 
INT18080 
INT18090 
INT18100 
INT18110 
INT18120 
INT18130 
INT18140 
INT18150 
INT18160 
INT18170 
INT18180 
INT18190 
INT18200 
INT18210 
INT18220 
INT18230 
INT18240 
INT18250 
INT18260 
INT18270 
INT18280 
INT18290 
INT183O0 
INT18310 
INT18320 
INT18330 
INT18340 


96 


20 
30 


10 
C 


c 
c 

c 


c 
c 
c 
c 
c 
c 


c 
c 

c 


JT  =  2 

DO  10  I  =  1,N2 

XM      =  X2(I) 

DO  20  J  =  JT,N1 

Jl  =  J  -1 

IF  (Xl(J).GE.XM 

CONTINUE 

J  =  Nl 

JT  =  J 

DXX     =  X2(I) 

DXX2    =  DXX  * 

DXX3    =  DXX2  * 

F2(I)    =  F1(J1) 


)  GO  TO  30 


•  X1(J1) 
DXX  /  2. 
DXX  /  3. 
+  DXX*FP1(J1) 


+  DXX2*FPP1(J1)  +  DXX3*FPPP1(J1) 


RETURN 
END 

DATA  SET 

DATA  SET 

DATA  SET 

SUBROUTINE  JO 

COMMON  /BLCO/ 

COMMON  /BLC1/ 

COMMON  /BLC7/ 

COMMON/EDDY 1/ 

I- 

COMMON  /GTY  / 
COMMON  /GRD  / 
COMMON  /SAVE/ 
COMMON  /SMRY/ 


AS  OF  08/24/84 
AS  OF  08/24/84 
AS  OF  02/20/84 


KCBCJOIN   AT  LEVEL  001 
KCBCJOIN   AT  LEVEL  001 
KCBCJOIN   AT  LEVEL  012 
IN( INDEX) 
NX,NXT,NP,NPT,NTR,IT,INVRS,NS,IP 
F(101,2),U(101,2),V(101,2)3W(101,2),B(101,2) 
C(100,100),D(100),DB(100),DBP'(100),UEO(100),GI 
RL , RX , S  QRX , RXNTR , GMTR , GMTR  S ( 1 0  0 ) 

,ALFAS( 100) ,FFS( 100) ,RTS( 100) , IEDY,NXSPT 
X(101),UE(100),P1(100),P2(100),CEL,CELH 
ETA( 101) ,DETA( 101) ,A( 101) 
FS(101),US(101),VS(101),BS(101),WS(101) 
VW( 100) , ITP( 100) , ISL( 100) ,DLS( 100) ,CF(  100) , 
THT(100),NPSTR(100) 


10 


15 


INDEX  =  1  FOR  THE  FIRST  SWEEP 
INDEX  =  2  FOR  SUBSEQUENT  SWEEP 

CALL  COMPGI 

CI I    =  C(NX,NX) 

UES    =  GI  /  (1.0  -  DLS(NX)  *  SQRT(RL)  *  CII  ) 

IF( INDEX. EQ. 1)  GOTO  15 

RETRIEVE  PROFILES  AT  STATION  NS  FOR  INVERSE  B.  L. 
CALCULATION 

DO  10  J=1,NPT 

F(J,2)  =  FS(J) 

U(J,2)  =  US(J) 

V(J,2)  =  VS(J) 

W(J,2)  =  WS(J) 

B(J,2)  =  BS(J) 

CONTINUE 

UES    =  UES/W(1,2) 

SQS    =  1.  0 

GOTO  30 

CONTINUE 

SQS    =  1. 0  /  SQRT(UES) 

DO  20  J=2,NPT 


INT18350 

INT18360 

INT18370 

INT1S380 

INT18390 

INT18400 

INT18410 

INT18420 

INT18430 

INT18440 

INT18450 

INT18460 

INT18470 

INT18480 

INT18490 

INT18500 

INT18510 

INT 185 20 

INT18530- 

INT18540 

INT18550 

INT18560 

INT18570 

INT18580 

INT18590 

INT18600 

INT18610 

INT18620 

INT18630 

INT18640 

INT18650 

INT18660 

INT18670 

INT18680 

INT18690 

INT18700 

INT18710 

INT18720 

INT18730 

INT18740 

INT18750 

INT18760 

INT18770 

INT18780 

INT18790 

INT18800 

INT18810 

INT18S20 

INT18830 

INT18840 

INT18850 

INT18860 

INT18870 

INT18880 

INT18890 

INT18900 


97 


20 

C 

30 


60 


C 

c 


65 

70 


80 


C 
C 

c 
c 
c- 


ETA(J)  =  ETA(  T")*SQS 
DETA(J-l)  =  ETA(J)  -  ETA(J- 
A(J)    =  0.5*DETA(J-1) 

CONTINUE 


1) 


DO  60  J=1,NPT 

U(J,2)  =  U(J,2)*UES 

UES  *  W(J,2) 

F(J,2)*SQS*UES 

V(J,2)/SQS*UES 


W(J,2)  = 

F(J,2)  = 

V(J,2)  = 

CONTINUE 

UE(NX)  = 

RX 

SQRX 

IF(NX. GT 

CALL  FILLUP(2) 

IF(  INDEX.  EQ.  2) 


W(l,2) 

UE(NX)*X(NX)*RL 

SQRT(RX) 

NTR)  CALL  EDDY 


GOTO  70 


STORE  PROFILES  AT  STATION  NS  FOR  THE  NEXT  SWEEP 
DO  65  J=1,NPT 
FS(J)   =  F(J,2) 
US(J)   =  U(J,2) 
VS(J)   =  V(J,2) 
WS(J)   =  W(J,2) 
BS(J)   =  B(J,2) 
CONTINUE 
DO  80  J=1,NPT 
F(J,1)  =  F(J,2) 
U(J,1)  =  U(J,2) 
V(J,1)  =  V(J,2) 
W(J,1)  =  W(J,2) 
B(J,1)  =  B(J,2) 
CONTINUE 
RETURN 
END 

DATA  SET  KCBCMAIN 

PROGRAM  MAIN 


AT  LEVEL  005  AS  OF  09/18/84 


SUBROUTINE  CASBLP( K2 , XP , YP , XMP , YMP , XS , YS , DLSP , VC , DBPP 
+  ,RN,NBL,ITRI,XCTRI,CASEID) 

COMMON  /BLCO/  NX,NXT,NP,NPT,NTR, IT, INVRS,NS, IP 

COMMON/BLIN/  TITLE(20) ,XC( 100) , YC( 100) , ISG( 100) ,DELS( 100) , 
+  XCTR,XTR,ISTRP,ICYCLE,ICYTL,XCTRS(2),TRFIND(2) 

C0MM0N/EDDY1/RL, RX, SQRX, RXNTR,GMTR,GMTRS( 100), 
+  ALFAS( 100) ,FFS( 100) ,RTS( 100) , IEDY,NXSPT 

COMMON/ BLOW/  VN(IOO) 

COMMON/TRN/  PGAMTR , OMEGA , RTHETB , RTRANB 

C0MM0N/PL0T/NVP(2),NXVP(20,2),ICC 


DIMENSION 

CASEID( 

20  ), 

XCTRI 

( 

2  ), 

ITRI 

( 

2  ) 

DIMENSION 

XP 

( 

100), 

YP 

( 

100), 

XMP 

( 

100) 

DIMENSION 

YMP 

( 

100), 

VC 

( 

100), 

SM 

( 

100) 

DIMENSION 

XS 

( 

100), 

YS 

( 

100), 

NBL 

( 

2  ) 

DIMENSION 

VT 

( 

100), 

S 

( 

100), 

DLSP 

( 

100) 

DIMENSION 

DLS 

( 

100), 

xo 

( 

100), 

YO 

( 

100) 

2 


INT18910 

INT18920 

INT18930 

INT18940 

INT18950 

INT18960 

INT18970 

INT18980 

INT18990 

INT19000 

INT19010 

INT19020 

INT19030 

INT19040 

INT19050 

INT1906f 

INT1907 

INT1908 

INT19090 

INT19100 

INT19110 

INT19120 

INT19130 

INT19140 

INT19150 

INT19160 

INT19170 

INT19180 

INT19190 

INT19200 

INT19210 

INT19220 

INT19230 

INT19240 

INT19250 

INT19260 

INT19270 

INT19280 

INT19290 

INT19300 

INT19310 

INT19320 

INT19330 

INT19340 

INT19350 

INT19360 

INT19370 

INT19380 

INT19390 

INT19400 

INT19410 

INT19420 

INT19430 

INT19440 

INT19450 

INT19460 


98 


DIMENSION 

DIMENSION 

c 

10 

FORMAT  ( 

20 

FORMAT  ( 

30 

FORMAT  ( 

40 

FORMAT  ( 

100 

FORMAT  ( 

110 

FORMAT  ( 

120 

FORMAT  ( 

+ 

130 

FORMAT  ( 

+ 

140 

FORMAT  ( 

150 

FORMAT  ( 

+ 

160 

FORMAT  ( 

170 

FORMAT  ( 

+ 

180 

FORMAT  (I 

C 

Dl     (  100),   D2     (  100),   D3     (  100) 
XPS    (100),   YPS(IOO)  ,  DBPP(IOO) 

20A4  ) 

415  ) 

3I5,3F10. 0  ) 

6F10.5  ) 

1H1,5X,20A4  ) 

1H0,6X, 'CYCLE  NO.  =  ' ,15  ) 

1H0,6X, 'FINISHED  TOTAL  NUMBER  OF  CYCLES  =  ',15, 

4X,'J0B  COMPLETED. '  ) 
lHO,*X,fTHE  UPDATED  DISPLACEMENT  SURFACE' ,/lX, 2HNX, 

9X , 2HXP , 12X , 2HYP , 1 IX , 3HDLS , 10X , 4HDBPP , 12X , 2HVN) 
I5,5E14.  6) 
1H0 . 6X , *  READ  IN  CONTROL  POINTS  DATA' , /1X, 2HNX, 9X, 

3HXMP,11X,3HYMP,12X,2HSM,13X,2HVC) 
I5,4E14.  6) 

1H0,6X, 'INTERPOLATED  ORIGINAL  GEOMETRY  DATA',/1X, 
2HNX , 9X , 2HXP , 12X , 2HYP , 12X , 1HS , 14X , 2HVT) 


C 
C 
5 


C 
C 
C 


15 


50 


60 
C 


GRANG(X1,X2,X3,Y1,Y2,Y3,X0)=  (X0-X2)*(X0-X3)/(X1-X2)/(X1-X3)*Y1 
+  +(X0-X1)*(X0-X3)/(X2-X1)/(X2-X3)*Y2+(X0-X1)*(X0-X2) 

+  /(X3-X1)/(X3-X2)*Y3 

ISWPT  =  1 

IEDY  =  1 
NBL(l)  =  91 
NBL(2)  =  71 

WRITE  (  6,100  )  CASEID 
CONTINUE 
ISTRP  =  I CYCLE 
NN  =  K2  -  1 
NT  =  K2 
WRITE  (  6,110  )  ICYCLE 
INTERPOLATE  OUTPUT  CONTROL  POINTS  TO  ORIGINAL  GEOMETRY 

DO  15  I  =  1  ,  NT 
XPS(I)  =  XP(I) 
YPS(I)  =  YP(I) 
CONTINUE 
S(l)  =  0.  0 
DO  50  I  =  2  ,  NT 

S(I)  =  S(I-l)  +  SQRT((XP(I)-XP(I-1))**2  + 
+  (YP(I)-YP(I-1))**2) 

CONTINUE 

SM(1)  =  SQRT((XMP(1)-XP(1))**2  +  ( YMP( 1) -YP( 1) )**2) 
DO  60  I  =  2  ,  NN 

SM(I)  =  SM(I-l)  +  SQRT((XMP(I)-XMP(I-1))**2+ 
+  (YMP(I)-YMP(I-1))**2) 

CONTINUE 

CALL  AMEAN(NN-10,NN,SM,VC,1) 
CALL  AMEAN(1,NN,SM,VC,1) 
SNT  =  S(NT) 


INT19470 
INT19480 
INT19490 
INT19500 
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INT19640 
INT19650 
INT19660 
INT19670 
INT19680 
INT19690 
INT19700 
INT19710 
INT19720 
INT19730 
INT19740 
INT19  750 
INT19760 
INT19770 
INT19780 
INT19790 
INT19800 
INT19810 
INT19820 
INT19830 
INT19840 
INT19850 
INT19860 
INT19870 
INT19880 
INT19890 
INT19900 
INT19910 
INT19920 
INT19930 
INT19940 
INT19950 
INT19960 
INT19970 
INT19980 
INT19990 
INT20000 
INT20010 
INT20020 


99 


65 


C 
C 
C 
C 
C 
C 
C 


44 


45 


C 
C 
C 
C 
C 
C 

c 
c 


70 
C 


80 

C 
C 


C 
C 
C 


SM(NT)  =  SM(NN)+SQRT( (XMP(NN) -XP(NT) )**2+(YMP(NN) -YP(NT) )**2) 

SMNT  =  S(NT)  /  SM(NT) 
DO  65  I  =  1  ,  NN 
SM(I)  =  SM(I)  *  SMNT 
CALL  DIFF3(NN,SM,VC,D1,D2,D3,0) 
CALL  INTRP3(NN, SM , VC ,D1 ,D2 ,D3 ,NT, S , VT) 
PRINT  OUT  INPUT  DATA 

WRITE(6,150) 

WRITE(6,160)  (I,XMP(I),YMP(I),SM(I),VC(I),I=1,NN) 

WRITE(6,170) 

WRITE(6,180)  (I,XP(I),YP(I),S(I),VT(I),I=1,NT) 

XPMIN  =  XP(1) 
DO  44  I  =  2  ,  NT 
IF  (XP(I)  . GT.  XPMIN)  GO  TO  44 
XPMIN  =  XP(I) 
CONTINUE 

CHORD  =  XP(NT)  -  XPMIN 
DO  45  I  =  1  ,  NT 
XP(I)  =  (XP(I)-XPMIN)  /  CHORD 
YP(I)  =  YP(I)  /  CHORD 
CONTINUE 

CALL  COMPBL  (  CASEID,XP, YP , VT, S ,DLSP,DLS ,DBPP,NBL, ITRI ,XCTRI , 
+  RN,NT,ISWPT) 


SMOOTH  THE  CALCULATED  DISPLACEMENT  THINKNESS 

CALL  SMFIT(1,NT,S,DLS,D1,2) 

ADD  DISPLACEMENT  THINKNESS  ON  THE  ORIGINAL  BODY 

DO  70  I  =  1  ,  NT 

DLSP(I)  =  DLS(I) 

CONTINUE 

CALL  SMFIT  ( 1 ,NT,XS , YS ,D1 , 2) 

DO  80  I  =  1  ,  NT 

XP(I)  =  XPS(I) 

YP(I)  =  YPS(I) 

CONTINUE 

IF  (ICYCLE  . GE.  ICYTL-1  .OR.  IP  . GE.  0)  THEN 

WRITE  (6,130) 

WRITE  (6,140)  (I,XP(I),YP(I),DLS(I),DBPP(I),VN(I),I=1,NT) 
WRITE  (  6,120  )  ICYCLE 
END  IF 

RETURN 
END 

DATA  SET  KCBCMAIN2  AT  LEVEL  001  AS  OF  08/24/84 
DATA  SET  KCBCMAIN2  AT  LEVEL  001  AS  OF  08/24/84 
DATA  SET  KCBCMAIN2  AT  LEVEL  010  AS  OF  04/06/84 
SUBROUTINE  MAIN2( ITR, ISWPT,SURFID) 
COMMON  /BLCO/  NX,NXT,NP,NPT,NTR, IT, INVRS ,NS , IP 
COMMON  /BLC1/  F( 101 ,2) ,U( 101 , 2) ,V( 101 ,2) ,W( 101 ,2) ,B( 101 ,2) 
COMMON  /BLC2/  DELF( 101) ,DELU( 101) ,DELV( 101) ,DELW( 101) 


INT20030 

INT20040 

INT20050 

INT20060 

INT20070J 

INT20080 

INT20090 

INT20100 

INT20110 

INT20120 

INT20130 

INT20140 

INT20150 

INT20160 

INT2017CJ 

INT20180 

INT20190 

INT20200 

INT20210 

INT20220 

INT20230 

INT20240 

INT20250 

INT20260 

INT20270I 

INT20280 

INT20290 

INT20300 

INT20310 

INT20320 

INT20330 

INT20340 

INT20350 

INT20360 

INT20370 

INT20380 

INT20390 

INT20400 

INT20410 

INT20420 

INT20430 

INT20440 

INT20450 

INT20460 

INT20470 

INT20480 

INT20490 

INT20500 

INT20510 

INT20520 

INT20530 

INT20540 

INT20550 

INT20560 

INT20570 

INT20580 


100 


c 
c 


10 


15 
20 

25 


30 


70 


COMMON  /BLC7/ 
COMMON  /BONV/ 
COMMON/ EDDY  1/ 


COMMON 
COMMON 


C(1005100),D(100),DB(100),DBP(100),UEO(100),GI 

ITMAX,EPSL,EPST,CONV 

RL,RX,SQRX,RXNTR,GMTR,GMTRS(100) 

,ALFAS( 100) ,FFS( 100) ,RTS( 100) , IEDY,NXSPT 
/GTY  /  X(101),UE(100),P1(100),P2(100),CEL,CELH 
/SMRY/  VW( 100) , ITP( 100) , ISL( 100) ,DLS( 100) ,CF( 100) , 
+  THT(100),NPSTR(100) 

COMMON  /BUN/  TITLE(  20)  ,XC(  100)  ,  YC(  100)  ,  ISG(  100)  ,DELS(  100)  , 
+  XCTR,XTR,ISTRP,ICYCLE,ICYTL,XCTRS(2),TRFIND(2) 

COMMON/BLC9/  UEB(IOO)  ,  CFS(IOO) 
COMMON/TRN/  PGAMTR , OMEGA , RTHETB , RTRANB 
COMMON  /I SURF/  ISF 
COMMON/PLOT/NVP( 2) ,NXVP( 20 , 2) , ICC 
DIMENSION  SURFID(4),RTSS(11) 
LOGICAL  SMOOTH  ,  SEPART  ,  HOMOPY 
DATA  RTSS/0.  0  ,0.  1 ,0.  2  ,0.  3,0.  4,0.  5  ,0.  6,0.  7  ,0.  8 ,0.  9  , 1.  0/ 


GRANG(X1,X2,X3,Y1,Y2,Y3,X0)=  (X0-X2)*(X0-X3)/(X1-X2)/(X1-X3)*Y1 


+ 

+ 


ISWP 

INDEX  = 
IGROWT  = 
NXSPT  = 


0 

1 
2 

NXT 


+CX0-X1  )*(  xo 
/(X3-X1)/(X3 


■X3)/(X2-X1)/(X2- 
■X2)*Y3 


■X3)*Y2+(X0-X1)*(X0-X2) 


+  1 


CALL  JOIN( INDEX) 


NXSTOP 

IF  (NS 

ISWP 

NX 

HOMOPY 

CEL 

Pl(NX) 

P2(NX) 

CELH 

IT 

CALL  COMPGI 

IGROW=l 

IT     =  IT  +  1 

RX     =  UE(NX)*X(NX) 

SQRX   =  SQRT(RX) 


GOTO  15 


=  NXT-1 

. GE.  NTR 

=  ISWP  +  1 

=  NX  +  1 

=  .FALSE. 

=  0.5*(X(NX)+X(NX- 

=  0.  5 

=  0.  0 

=  0.  5*CEL 

=  0 


1))/(X(NX)-X(NX-1)) 


'RL 


IF(IT  .  LE.  ITMAX)  GO  TO  80 

1F( HOMOPY)  GO  TO  72 

IRC  =  1 

RT  =  RTSS(IRC) 

HOMOPY  =  .TRUE. 

UEREF  =  UEO(NX-l) 

UESAVE  =  UEO(NX) 

UEO(NX)  =  RT*UESAVE+(1. 0-RT)*UEREF 

DO  61  J=1,NP 

F(J,2)  =  F(J,1) 

U(J,2)  =  U(J,1) 

V(J,2)  =  V(J,1) 

W(J,2)  =  W(J,1) 


INT20590 
INT20600 
INT20610 
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INT20640 
INT20650 
INT20660 
INT20670 
INT20680 
INT20690 
INT20700 
INT20710 
INT20720 
INT20730 
INT20740 
INT20750 
■INT20760 
INT207  70 
INT20780 
INT20790 
INT20800 
INT20810 
INT20820 
INT20830 
INT20840 
INT20850 
INT20860 
INT20870 
INT20880 
INT20890 
INT20900 
INT20910 
INT20920 
INT20930 
INT20940 
INT20950 
INT20960 
INT20970 
INT20980 
INT20990 
INT21000 
INT21010 
INT21020 
INT21030 
INT21040 
INT21050 
INT21060 
INT21070 
INT21080 
INT21090 
INT21100 
INT21110 
INT21120 
INT21130 
INT21140 


61 

C 
72 

C 
C 


80 

C 
C 
C 


C 
C 
C 
C 


B(J,2)  =  B(J,1) 

CONTINUE 

GO  TO  30 

NXSTOP  =  NX  -  1 

CALL  AMEAN( NS , NXSTOP , X , CF , 1 ) 

CALL  AMEAN( NS, NXSTOP, X,VW,1) 
CALL  HEADER(  TITLE , SURFID, ISTRP  ) 
WRITE  (6,  250  )  ISWP 

WRITE  (6,  260  )  (M,XC(M) ,X(M) ,CF(M) ,DLS(M) ,THT(M) ,UE(M) , 
+    UEO(M) ,D(M) ,DB(M) ,GMTRS(M) , ITP(M) ,NPSTR(M) ,M=1, NXSTOP) 
WRITE(6,  270  )  NX 
STOP 

CONTINUE 
IF(NX  .GT.  NTR)  GOTO  100 

LAMINAR  FLOW  CALCULATION 

CALL  C0EF(GAMMA1,GAMMA2) 

CALL  S0LV4 ( GAMMA 1 , GAMMA2 ) 

UE(NX)  =  U(NP,2) 

IF(ABS(DELV(1))  . GT.  EPSL)  GO  TO  70 

CHECK  ON  LAMINAR  FLOW  SEPARATION.  IF  SEPARATION  OCCURS,  ASSIGN  BEGIN 
OF  TRANSITION  TO  THAT  POINT  AND  RECOMPUTE  THE  CURRENT  STATION  NX 


IF(V(1,2).GT.  0.  0 
CALL  TRNS(ICODE) 
GOTO  25 


OR.  ITR.NE.  3)  GOTO  110 


C 
C 
C 
100 


TURBULENT  FLOW  CALCULATION 

CONTINUE 

CALL  EDDY 

CALL  COEF( GAMMA 1,GAMMA2) 

CALL  S0LV4( GAMMA 1,GAMMA2) 

UE(NX)  =  U(NP,2) 

VM     =  AMAX1(V(1,2),1.0) 

IF(ABS(DELV(1)/VM)  .  GT.  EPST)  GO  TO  70 


GROWTH 


110 

CONTINUE 

c 

c 

CHECK  FOR  B.  L. 

c 

IF(NP  .GE.  NP7) 

IF(ABS(V(NP,2)) 

+ 

CALL  FILLUPO) 

IGROW=IGROW+l 

IT   •   =  1 

GO  TO  70 

c 

120 

CONTINUE 

CALL  FILLUP(2) 

CALL  OUTPUTC2) 

IF(NX.  GE.NTR   . 

GO  TO  120 

.LT.  0.0005  .AND.  ABS(  1.  0-U(NP-2,2)/U(NP,2)) 

LT.  0.0035.  OR.  IGROW.GT.  IGROWT)  GOTO  120 


OR.  ITR.  EQ.  0)  GOTO  150 
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c 
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c 
c 
c 
c 


IF(NX.  LT.  3  .OR.   ITR.NE.  3)  GOTO  150 

CALCULATE  TRANSITION  LOCATION  USING  MICHEL  METHOD 

CALL  TRNS(ICODE) 
IF(ICODE.EQ.  0)  GOTO  150 

TRANSITION  OCCURS  BASED  ON  MICHEL  CRITERIOR  AT  STATION  NX 
RECALCULATE  B.  L.  AT  NX  STATION  ASSUMING  THE  FLOW  IS  TRANSITIONAL 


IT      =0 
I GROW   =  1 
GOTO  70 
150   CONTINUE 

IF(.NOT.  HOMOPY  )  GO  TO  154 
IF(  RT  .  GT.  0.9999)  GO  TO  154 
IRC  =  IRC  +  1 
RT  =  RTSS(IRC) 

UEO(NX)  =RT*UESAVE  +  ( 1. 0-RT)*UEREF 
GO  TO  30 
154   CONTINUE 

IF(NX  . LT.  NXSTOP)  GO  TO  20 
C 

C   THE  B.  L.  CALCULATION  FOR  THE  CURRENT  SWEEP  IS  COMPLETED. 
C   CHECK  FOR  THE  CONVERGENCE  AND  ,  IT  NOT,  MOVE  TO  THE  NEXT 
C    SWEEP. 
C 
C 
160   CONTINUE 

D(NXT)   =  GRANG(X(NXT-3) ,X(NXT-2) ,X(NXT-1) ,D(NXT-3) ,D(NXT-2) , 
+  D(NXT-1),X(NXT)) 

DLS(NXT)=  GRANG(X(NXT-3) ,X(NXT-2) ,X(NXT-1) ,DLS(NXT-3) ,DLS(NXT-2) , 
+  DLS(NXT-1),X(NXT)) 

UE(NXT)  =  GRANG(X(NXT-3),X(NXT-2),X(NXT-l),UE(NXT-3), 
+  UE(NXT-2),UE(NXT-1),X(NXT)) 

DO  165  I  =  1  ,  NXSTOP 
165   CFS(I)  =  CF(I) 

CALL  AMEAN( NS, NXSTOP, X,CF,1) 
C  CALL  AME AN(NS, NXSTOP, X,VW,1) 
C     CALL  HEADER(  TITLE ,SURFID, ISTRP  ) 

IF(ICYCLE  -LT.  ICYTL-1  .AND.  IP  .  LT.  0)GO  TO  170 
WRITE  (6,  250  )  ISWP 

WRITE  (6,  262  )  1 ,XC( 1) ,X( 1) ,CF( 1) ,DLS( 1) ,THT( 1) ,UE( 1) , 
+  UE0(1),0. 0,GMTRS(1),ITP(1),NPSTR(1) 

WRITE  (6,  264  )  (M,XC(M) ,X(M) ,CF(M) ,DLS(M) ,THT(M) ,UE(M) , 
+     UE0(M),DLS(M)/THT(M),GMTRS(M),ITP(M),NPSTR(M),M=2, NXSTOP) 
IF  ((ICYCLE.EQ.  ICYTL).  AND.  (IP.EQ.  -2).  AND.  (NVP(ISF).NE.  0))  THEN 
WRITE( 12,800)  NS+1,NTR,XCTR 

WRITEC 12,810)  (XC(M),X(M),UE(M),CF(M),GMTRS(M),ISG(M), 
2  M=l, NXSTOP) 

800     FORMAT(2I5,F10.  6) 
810     FORMAT(5E15.5,I5) 
END  IF 
170    CONTINUE 
C 

DMAX   =  D(l) 
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180 
C 

c 
c 
c 


190 
195 

200 

205 


210 


250 


260 


262 


264 
270 


C 
C 

C 


DDMAX  =  ABS(D(1)  -  DB(1)) 

DO  ISO  I  =  2,NXT 

DMAX   =  AMAX1(  DMAX,D(I)  ) 

DD     =  ABS(D(I)  -  DB(I)) 

DDMAX  =  AMAX1(  DDMAX, DD  ) 

CONTINUE 

IF  (  ABS(  DDMAX  /  DMAX  )  . LE.  0.0050  )  RETURN 


UPDATE  D  FOR  THE  NEXT  SWEEP 

IF  (ISWP  .  GT.  1)  GO  TO  195 

DO  190  I  =  NS  ,  NXT 

D(I)  =  D(I)*(1.  0-K)MEGA*(UE(I)/UEO(I)-l. 0)) 

GO  TO  205 

IF  (ISWP  .EQ.  2)  GOTO  205 

DO  200  I  =  NS  ,  NXT 

D(I)  =  D(I)  *  (1. 0+OMEGA*(UE(I)/UEB(I)-1.  0)) 

IFflSWP  .  GE.  ISWPT  )  RETURN 

NX      =  NS 

NP      =  NPSTR(NX) 

INDEX   =  2 

DO  210  1=  1,NXT 

DB(I)    =  D(I) 

UEB(I)  =  UE(I) 

CONTINUE 

GOTO  10 

F0RMAT(1H0,'  **  SUMMARY  OF  INVERSE  BOUNDARY  LAYER  SOLUTIONS.  **'  , 
+  1H0,4X,'ISWP  =',I3/) 

FORMAT( 1H0 , 4X , 2HNX , 5X , 3HX/C , 9X , 1HX , 9X , 2HCF ,  8X , 
+  3HDLS,8X,3HTHT,9X,2HUE,  8X, 3HUE0 , 10X, 1HD  ,9X,2HDB,3X, 

+       4HGMTR,4X,2HIT,1X,2HNP/(1H  , 3X, 13 ,F10.  5 , 8E11.  4 ,F8.  4, 213) ) 

F0RMAT(1H0,4X,2HNX,6X,3HX/C,11X,1HX,10X,2HCF,9X, 
+  3HDLS , 9X , 3HTHT , 10X , 2HUE , 9X , 3HUEO , 1 IX , 1HH , 8X , 

+       4HGMTR,4X,2HIT,1X,2HNP/(1H  ,3X, 13, 9E12.  4,213) ) 

F0RMAT(1H  ,3X,I3,9E12. 4,213) 

F0RMAT(1H0,'  **  ITERATIONS  EXCEEDED  ITMAX  AT  NX  =',15,' 


.<.  i 


END 


1H0,'  **  CALCULATIONS  STOP.  **') 


,/ 


OF  08/24/84 
OF  08/24/84 
OF  02/22/84 


DATA  SET  KCBCOUTPUT  AT  LEVEL  001  AS 
DATA  SET  KCBCOUTPUT  AT  LEVEL  001  AS 
DATA  SET  KCBCOUTPUT  AT  LEVEL  002  AS 

SUBROUTINE  OUTPUT( INDEX) 

COMMON  /BLCO/  NX,NXT,NP,NPT,NTR, IT, INVRS ,NS , IP 

COMMON/BLIN/  TITLE(20) ,XC( 100) ,YC( 100) , ISG( 100) ,DELS( 100) , 
+  XCTR,XTR,ISTRP,ICYCLE,ICYTL,XCTRS(2),TRFIND(2) 

COMMON  /BLC1/  F( 101 , 2) ,U( 101 ,2) , V( 101 , 2) ,W( 101 ,2) ,B( 101 , 2) 

COMMON  /BLC7/  C( 100 , 100) ,D( 100) ,DB( 100) ,DBP( 100) ,UE0( 100) ,GI 

COMMON/EDDY1/  RL,RX,SQRX,RXNTR,GMTR,GMTRS( 100) 
+  ,ALFAS( 100) ,FFS( 100) ,RTS( 100) , IEDY,NXSPT 

COMMON  /GTY  /  X( 101) ,UE( 100) ,P1( 100) ,P2( 100) ,CEL,CELH 

COMMON  /GRD  /  ETA( 101) ,DETA( 101) , A( 101) 

COMMON  /SMRY/  VW( 100) , ITP( 100) , ISL( 100) ,DLS( 100) ,CF( 100) ,THT( 100) 
+  NPSTR(IOO) 

COMMON  /I SURF/  ISF 
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C0MM0N/PL0T/NVP( 2) ,NXVP( 20 , 2) , ICC 


5 
C 
C 
10 


20 


C 
C 

100 


120 


150 


ITP(NX)  =  IT 

NPSTR(NX)=NP 

IF(NX. GT. 1)  GOTO  5 

DLS(NX)=  0.  0 

VW(NX)  =  0.  0 

D(NX)   =  0.  0 

THT(NX)=  0. 0 

CF(NX)  =  0.  0 

VW(NX)  =  0.  0 

GOTO  150 

GOTO  (10,100,200),  INDEX 

CALCULATE  B.  L.  PARAMETERS  FOR  TRANSFORMED  COORDINATES 
CONTINUE 

CF(NX)  =2.0  *  V(l,2)  *  B(1,2)/SQRX 
VW(NX)  =  UE(NX)*SQRT(UE(NX)/X(NX))*V(1,2) 
DLS(NX)=  X(NX)/SQRX  *  (ETA(NP) -F(NP, 2) ) 
D(NX)   =  UE(NX)  *  DLS(NX)  *  SQRT(RL) 
Ul     =  U(l,2)  *  (1.0  -U(l,2)) 
SUM    =0.0 
DO  20  J=2,NP 

U2     =  U(J,2)  *  (1.  0  -U(J,2)) 
SUM    =  SUM  +  A( J)  *  (Ul  +  U2) 
Ul     =  U2 
CONTINUE 

THT(NX)=  X(NX)/SQRX  *  SUM 
GOTO  150 

CALCULATE  B.  L.  PARAMETERS  FOR  SEMI-TRANSF  COORDINATES 
CONTINUE 

SQXC  =  SQRT(X(NX)) 
SQRL  =  SQRT(RL) 

CF(NX)  =2.0  *  V(l,2)  *  B(1,2)/(SQXC*SQRL*W(NP,2)**2) 
VW(NX)  =  V(l,2)  /  SQXC 
UE(NX)  =  U(NP,2) 
RX  =  RL  *  UE(NX)  *  X(NX) 

DLS(NX)  =  (ETA(NP)-F(NP,2)/U(NP,2))/SQRL*SQXC 
SUM     =0.0 

Ul      =  U(1,2)/U(NP,2)*(1.0  -U(1,2)/U(NP,2)) 
DO  120  J=2,NP 

U2      =  U(J,2)/U(NP, 2)*(1.  0  -U(J,2)/U(NP,2)) 
SUM     =  SUM  +  A(J)  *  (Ul  +  U2) 
Ul      =  U2 
CONTINUE 

THT(NX)  =  SUM  /SQRL  *  SQXC 

D(NX)    =  (U(NP,2)--'"ETA(NP)-F(NP,2))  *  SQXC 
IF  (NX  .GE.  NXT)  GO  TO  160 
IF  (IEDY  .EQ.  0  .OR.  NX  .  LE.  NTR+2)  GO  TO  160 


C  MODIFY  ALFA  USING  SIMPSON'S  ARGUMENTS 
C 
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160 


175 

C 
C 

200 


C 
210 


4001 
4000 

4100 
4200 
4300 

C 
C 
C 

C 
C 
C 
C 
C 
C 
C 


CALL  SMPSON 

DO  175  J=1,NPT 

F(J,1)   =  F(J,2) 

U(J,1)   =  U(J  2) 

V(J,1)   =  V(J,2) 

W(J,1)   =  W(J,2) 

B(J,1)   =  B(J.2^ 

CONTINUE 

IF  ((IP.LE.O).  AND.  ((IP.NE. -2).0R.  (ICYCLE.LT.  ICYTL)))  RETURN 

PRINT  OUT  VELOCITY  PROFILES 
IF  (NX.  EQ.  1)   GOTO  210 
IF  (NX.  LE.NS)  THEN 

FAC1  =  SQRT(X(NX)/RL/UE(NX)) 

FAC2  =1.0 
ELSE 

FAC1  =  SQRT(X(NX)/RL) 

FAC2  =  1. 0/UE(NX) 
ENDIF 

NPM1   =  NP  -1 
WRITE(6,4001)  NX,X(NX) 
WRITE(6,4000) 

WRITE(6,4100)  (J,ETA(J),F(J,2),U(J,2),V(J,2),W(J,2),B(J,2), 
+  ETA(J)*FAC1,U(J,2)*FAC2,J=1,NPM1,3) 

WRITE(6,4100)  NP,ETA(NP),F(NP,2),U(NP,2),V(NP,2),W(NP,2),B(NP,2), 
+  ETA(NP)*FAC1,U(NP,2)*FAC2 

IF  (IP.NE.  -2)  RETURN 

IF  ((NXVP(ICC,ISF).NE.NX).OR.  (ICC.GT.  NVP(ISF)))  RETURN 

WRITE( 12,4200)  NP 

WRITE( 12,4300)  (ETA( J) , J=l ,NP) 

WRITE(12,4300)  (U( J, 2) , J=1,NP) 

ICC  =  ICC+1 

RETURN 

FORMAT(/1HO,'NX  =',15,'   S/C  =',F10.5) 

FORMAT( 1H0 ,2H  J,9X, 3HETA, 15X, 1HF, 13X, 1HU, 13X, 1HV, 13X, 1HV, 13X, 1HB , 
+       13X,3HY/C,10X,4HU/UE) 

F0RMATC1H  ,I3,E14. 5,2X,5E14. 5,2X,2E14.  5) 

F0RMAT(I5) 

FORMAT(8F10.6) 

END 

DATA  SET  KCBCSMFIT  AT  LEVEL  001  AS  OF  08/24/84 
DATA  SET  KCBCSMFIT  AT  LEVEL  001  AS  OF  08/24/84 
DATA  SET  KCBCSMFIT  AT  LEVEL  001  AS  OF  08/15/83 

SUBROUTINE  SMFIT(NS,ND,X,Q,D,KS) 

THIS  SUBROUTINE  SMOOTHES  DATA,  Q,  USING  FIVE -POINT  FORMULA. 


NS   :  BEGINNING  POINT?   ND 
X   :  INDEDENPENT  COORDINATE?  D 
Q   :  VARIABLE  TO  BE  SMOOTHED 
KS   :  NO  OF  SMOOTHING 

DIMENSION  X(101),Q(101),D(101) 


END  POINT 

WORKING  STORAGE 


SMT5(Q1,Q2,Q3,Q4,Q5    )  =  0. 0625*( 10. 0*Q3+4. 0*(Q2+Q4) -Q1-Q5) 
SMT3(Q1,Q2,Q3,X1,X2,X3)  =  0. 5*(Q2+(Q1»ABS(X3-X2)+Q3*ABS(X2-X1) ) 
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20 

40 
100 

C 
200 


220 


250 
300 
C 


C 
C 
C 


/ABS(X3-X1)) 

IF(KS.LE.  0)  RETURN 
NSP1    =  NS+1 
NSP2   =  NS+2 
NDM1   =  ND-1 
NDM2   =  ND-2 

NDIF   =  ND-NS+1 

IF  (  NDIF  . LT.  3  )  RETURN 

IF  (  NDIF  . LT.  5  )  GO  TO  200 

DO  100  K=1,KS 

D(NS+1)=  SMT3fQ(NS) ,Q(NS+1) ,Q(NS+2) ,X(NS) ,X(NS+1) ,X(NS+2)) 

D(ND-1)=  SMT3'  0(ND-2)JQ(ND-1),Q(ND)SX(ND-2),X(ND-1)JX(ND)) 

DO  20  I=NSP2,NDM2 

D(I)   =  SMT5(Q(I-2),Q(I-l),Q(I),Q(I+l),Q(I+2)) 

CONTINUE 

DO  40  I=N£ri,NDMl 

Q(I)     =  D(I) 

CONTINUE 

RETURN 

DO  300  K  =  1,KS 

DO  220  I  =  NSP1,NDM1 

D(I)    =  SMT3(Q(I-1),Q(I)JQ(I+1),X(I-1),X(I)JX(I+1)) 

CONTINUE 

DO  250  I  =  NSP1,NDM1 

Q(D    =  D(I) 

CONTINUE 

CONTINUE 


RETURN 
END 

DATA 


DATA 


AT  LEVEL 

001 

AS 

OF 

08/24/84 

AT  LEVEL 

001 

AS 

OF 

08/24/84 

AT  LEVEL 

005 

AS 

OF 

02/21/84 

+ 


+ 


SET  KCBCSOLV3 
SET  KCBCSOLV3 
DATA  SET  KCBCS0LV3 
SUBROUTINE  SOLV3 
COMMON  /BLCO/  NX,NXT,NP,NPT,NTR, IT, INVRS ,NS , IP 

F(101,2),U(101,2),V(101,2),W(101,2),B(101,2) 
DELF( 101) ,DELU( 101) ,DELV( 101) ,DELW( 101) 
S1(101),S2(101),S3(101),S4(101),S5(101),S6(101), 
S7(101),S8(101),R1(101),R2(101),R3(101),R4( 
ETA(101),DETA(101),A(101) 
A11(101),A12(101),A13(101),A14(101), 

A21(101),A22(101),A23(101),A24(101) 


COMMON 
COMMON 
COMMON 

COMMON 
COMMON 


/BLC1/ 
/BLC2/ 
/BLC6/ 

/GRD  / 
/BLCB/ 


All(l)=  1.0 
A12(l)=  0.0 
A13(l)=  0.  0 
A21(l)=  0.  0 
A22(l)=  1.  0 
A23(l)=  0.  0 
Gll=-1.  0 
G12=-A(2) 
G13=  0.  0 
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c 
c 
c 


G21=  S4(2) 
G23=-S2(2)/A(2) 
G22=  G23+S6(2) 

FORWARD  SWEEP 

DO  500  J=2,NP 


IF(J  . EQ.  2)  GO  TO  100 

DEN   =  (A13(J-1)*A21(J-1)-A23(J-1)*A11(J-1)-A(J)* 
+  (A12(J-1)*A21(J-1)-A22(J-1)*A11(J-1))) 

DEN1  =  A22(J-1)*A(J)-A23(J-1) 

Gll=  (A23(J-1)+A(J)*(A(J)*A21(J-1)-A22(J-1)))/DEN 

G12=  -(A(J)*A(J)+G11-'KA12(J-1)*A(J)-A13(J-1)))/DEN1 

G13=  ( Gl  1*A13(  J- 1  )+G12*A23(  J- 1 )  )  /A(  J) 

G21=  (S2(J)*A21(J-1)-S4(J)*A23(J-1)+A(J)*(S4(J)* 
+  A22(J-1)-S6(J)*A21(J-1)))/DEN 

G22=  (-S2(J)+S6(J)*A(J)-G21*(A(J)*A12(J-1)-A13(J-1)))/DEN1 

G23=  G21*A12(J-1)+G22*A22(J-1)-S6(J) 
100    A11(J)=  1.  0 

A12(J)=-A(J)-G13 

A13(J)=  A(J)*G13 

A21(J)=  S3(J) 

A22(J)=  S5(J)-G23 

A23(J)=  S1(J)+A(J)*G23 

R1(J)  =  Rl(J)-(Gll*Rl(J-l)+ul2*R2(J- 

R2(J)  =  R2(J)-'G21*R1(J-1)+G22*R2(J- 

CONTINUE 


1)+G13*R3(J-1)) 
<1)+G23*R3(J-1)) 


500 

C 

C   BACKWARD  SWEEP 

C 

DELU''NP)  =  R3(NP) 

El       =  R1(NP)-A12(NP)*DELU(NP) 

E2       =  R2(NP)-A22(NP)*DELU(NP) 

DELV(NP)  =  (E2*A11(NP)-E1*A21(NP))/(A23(NP)*A11(NP)-A13(NP) 
+  A21(NP)) 

DELF(NP)  =  (E1-A13(NP)--DELV(NP))/A11(NP) 

J     =  NP 
600   J     =  J-l 

E3    =  R3(J)-DELU(J+1)+A(J+1)*DELV(J+1) 

DEN2  =  A21(J)^A12(J)^A(J+1)-A21(J)->A13(J)-A(J+1)"A22(J)^A1 
+  A23(J)*A11(J) 

DELV(J)   =  (A11(J)*(R2(J)+E3*A22(J))-A21(J)*R1(J)-E3*A21(J) 
+  )/DEN2 

DELU(J)   =-A(J+l)*DELV(J)-E3 

DELF(J)   =  (R1(J)-A12(J)*DELU(J)-A13(J)*DELV(J))/A11(J) 

IF(J  .GT.  1)  GO  TO  600 
C 

DO  700  J=1,NP 

F(J,2)=  F(J,2)+DELF(J) 

U(J,2)=  U(J,2)+DELU(J) 

V(J,2)=  V(J,2)+DELV(J) 
700   CONTINUE 

U(l,2)=  0.  0 

CALL  EDGCHK(NP,ETA,F(1,2),U(1,2),V(1,2)) 

RETURN 
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END 
C  DATA  SET  KCBCS0LV4   AT  LEVEL  001  AS  OF  08/24/84 

C  DATA  SET  KCBCS0LV4  AT  LEVEL  001  AS  OF  08/24/84 

C  DATA  SET  KCBCS0LV4   AT  LEVEL  001  AS  OF  02/21/84 

SUBROUTINE  S0LV4( GAMMA 1 .GAMMA2) 
COMMON  /BLCO/  NX,NXT,NP,NPT,NTR, IT, INVRS ,NS , IP 
COMMON  /BLC1/  F( 101 , 2) ,U( 101 , 2) , V( 101 , 2) ,W( 101 ,2) ,B( 101 ,2) 
COMMON  /BLC2/  DELF( 101) ,DELU( 101) ,DELV( 101) ,DELW( 101) 
COMMON  /BLC6/  Sl( 101) , S2( 101) , S3( 101) , S4( 101) , S5( 101) , S6( 101)  , 
+  S7(101),S8(101),R1(101),R2(101),R3(101),R4( 

COMMON  /GRD  /  ETA( 101) ,DETA( 101) ,A( 101) 
COMMON  /BLCB/  All( 101) , A12( 101) ,A13( 101) ,A14( 101) , 
+  A21(101),A22(101),A23(101),A24(101) 


All(l) 

=  1.  0 

A12(l) 

=  0.0 

A13(l) 

=  0.0 

A14(l) 

=  0.0 

A21(l) 

=  0.0 

A22(l) 

=  1.  0 

A23(l) 

=  0.0 

A24(l) 

=  0.  0 

DO  10  J 

=  2,NP 

AA1 

=  A13(J-1)-A(J)*A12(J-1) 

AA2 

=  A23(J-1)-A(J)*A22(J-1) 

AA3 

=  S2(J)-A(J)*S6(J) 

DET 

=  AA2*A11(J-1)-AA1*A21(J-1) 

AJS 

=  A(J)**2 

Gil 

-(AA2+A21(J-1)*AJS)/DET 

G12 

(A11(J-1)*AJS+AA1)/DET 

G13 

A12(J-1)*G11+A22(J-1)*G12+A(J) 

G14 

A14(J-1)*G11+A24(J-1)*G12 

G21 

(S4(J)*AA2-A21(J-1)*AA3)/DET 

G22 

( Al  1  ( J- 1  )*AA3  -S4(  J)*AA1 )  /DET 

G23 

A12(J  •.)*G21+A22(J-1)*G22-S6(J) 

G24 

A14(J-1)*G21+A24(J-1)*G22-S8(J) 

All(J) 

=  1.0 

A12(J) 

=  -A/  -)-G13 

A13(J) 

=  A(J)^G13 

A14(J) 

=  -G14 

A21(J) 

=  S3(J) 

A22(J) 

=  S5(J)-G23 

A23(J) 

=  S1(J)+A(J)*G23 

A24(J) 

=  S7(J)-G24 

R1(J) 

=  R1(J)  -G11*R1(J-1)-G12*R2(J- 

+ 

-G14-'-R4(J-l) 

R2(J) 

=  R2(J)  -G21*R1(J-1)-G22*R2(J- 

+ 

-G24*R4(J-1) 

10    CONTINUE 

r 

C   BACKWARD 

SWEEP 

J 

=  NP 

Gl 

=  GAMMA 1/GAMMA2 

R3(J) 

=  R3(J)/GAMMA2 

1)-R3(J-1)*G13 
1)-R3(J-1)*G23 
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Rl(J) 

R2(J) 

CI 

C2 

DET 

DELF(J)   : 

DELV(J)   = 

DELW(J)  ■■ 

DELU(J)   = 
20    J 

CC1 

CC2 

CC3 

CC4 

CC5 

CC6 

DENO 

DELF(J)   = 

DELV(J)   = 

DELV(J)   = 

DELU(J)   : 

IF(J  .GE. 

DO  30  J  = 

F(J,2)  = 

U(J,2)  = 

V(J,2)  = 

W(J,2)  = 
30    CONTINUE 

U(l,2)  =  0.  0 

CALL  EDGCHK(NP,ETA,F(1,2),U(1,2),V(1,2)) 

RETURN 


C 
C 

c 

c 
c 
c 


c 

100 
110 

c 


=  R1(J)-A12(J)*(R4(J)+R3(J))-A14(J)*R3(J) 

=  R2(J)-A22(J)*(R4(J)+R3(J))-A24(J)*R3(J) 

=  A11(J)-G1*(A12(J)+A14(J)) 

=  A21(J)-G1*(A22(J)+A24(J)) 

=  C1*A23(J)-C2*A13(J) 

=  (R1(J)*A23(J)-R2(J)*A13(J))/DET 

=  (C1*R2(J)-C2*R1(J))/DET 

=  R3(J)-G1*DELF(J) 

=  R4(J)+DELW(J) 

=  J-l 

=  DELU(J+1)-R3(J)-A(J+1)*DELV(J+1) 

=  DELW(J+1)-R4(J) 

=  A13(J)-A(J+1)*A12(J) 

=  R1(J)-A12(J)*CC1-A14(J)*CC2 

=  A23(J)-A(J+1)*A22(J) 

=  R2(J)-A22(J)*CC1-A24(J)*CC2 

=  A11(J)*CC5-A21(J)*CC3 

=  (CC4*CC5-CC3*CC6)/DENO 

=  (All(J)*CC6-A21(J)*CC4)/DENO 

=  CC2 

=  CC1-A(J+1)*DELV(J) 

2)  GO  TO  20 

1,NP 
F(J,2)+DELF(J) 
U(J,2)+DELU(J) 
V(J,2)+DELV(J) 
W(J,2)+DELW(J) 


END 

DATA  SET  KCBCTRNS 
DATA  SET  KCBCTRNS 
DATA  SET  KCBCTRNS 

SUBROUTINE  TRNS(ICODE) 


AT  LEVEL  001  AS  OF  08/24/84 
AT  LEVEL  001  AS  OF  08/24/84 
AT  LEVEL  005  AS  OF  03/13/84 


CALCULATE  TRANSITION  LOCATION  USING  MICHEL  CRITERION 

COMMON  /BLCO/  NX,NXT,NP ,NPT,NTR, IT, INVRS ,NS , IP 

COMMON  /BLC1/  F( 101,2) ,U( 101,2) ,V( 101 ,2) ,W( 101,2) ,B( 101,2) 

COMMON/BLIN/  TITLE( 20) ,XC( 100) ,YC( 100) , ISG( 100) ,DELS( 100) , 

+  XCTR,XTR,ISTRP,ICYCLE,ICYTL,XCTRS(2),TRFIND(2) 

COMMON/ EDDY 1/  RL , RX , SQRX , RXNTR , GMTR , GMTRS( 100 ) 

+  ,ALFAS( 100) ,FFS( 100) ,RTS( 100) ,IEDY,NXSPT 

COMMON  /GRD  /  ETA( 101) ,DETA( 101) , A( 101) 
COMMON  /GTY  /  X( 101) ,UE( 100) ,P1( 100) ,P2( 100) ,CEL,CELH 
COMMON/TRN/  PGAMTR , OMEGA , RTHETB , RTRANB 

FORMAT(/3X. '2EG1N  OF  TRANSITION  HAS  BEEN  DETECTED  BY  MICHEL"  S  ' 
+' CRITERION:  X/C  =' ,F8. 4,4X, ' S/C  =' ,F8. 4,4X, 'NTR  =',I3/) 

F0RMAT(/3X,' BEGIN  OF  TRANSITION  IS  ASSUMED  AT  THE  POINT  OF  * , 
+' LAMINAR  SEPARATION:  X/C  =' ,F8. 4,4X, ' S/C  =' ,F8. 4,4X, ' NTR  =',I3/) 


I  CODE 


=  0 
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ISEP 


=  1 


C 

c 


c 
c 


c 
c 

c 

c 


IF(V(1,2).  LT.  0.  0)  THEN 
***  TRANSITION  PROCESS  HAS  BEGUN  DUE  TO  LAMINAR  SEPARATION  *** 
FAC  =  V(1,1)/(V(1,1)-V(1,2)) 
GOTO  20 
END  IF 

***  CHECK  MICHEL'S  TRANSITION  CRITERION  *** 
ISEP   =  0 
SUM    =  0.  0 

Fl     =  U(1,2)/U(NP,2)*(1.0-U(1,2)/U(NP,2)) 
DO  10  J=2,NP 

F2     =  U(J,2)/U(NP,2)*(1.0-U(J,2)/U(NP,2)) 
SUM    =  SUM  +  (Fl  +  F2  )  *A(J) 
Fl     =  F2 
10  CONTINUE 

CONV   =  SQRT(RL/X(NX)) 

IF(NX. LE.NS)  CONV  =  SQRT(RX)/X(NX) 

THETA  =  SUM  /  CONV 

RTHETA  =  RL  *  UE(NX)  *  THETA 

RTRAN  =  1.174  *  ( 1.  0+22400. 0/RX)  *  RX**0. 46 

IF( RTHETA.  LT. RTRAN)  THEN 

RTHETB  =  RTHETA 

RTRANB  =  RTRAN 

RETURN 
END  IF 

***  TRANSITION  PROCESS  HAS  BEGUN  BECAUSE  OF  MICHEL'S  CRITERION  *** 
FAC  =  (RTHETB-RTRANB)/(RTRAN-RTRANB-RTHETA+RTHETB) 

***  COMPUTE  EXACT  LOCATION  OF  TRANSITION  BEGIN  *** 
20  NTR  =  NX-1 
NTR1  =  NTR  +  1 

XCTR  =  XC(NX-l)  +  FAC*(XC(NX)-XC(NX-1)) 
XTR  =  X(NX-l)  +  FAC*(X(NX)-X(NX-1)) 
UETR  =  UE(NX-l)  +  FAC--'-(UE(NX)-UE(NX-l)) 
IF  (  ISEP  .EQ.  0  )  WRITE  (6,100)  XCTR, XTR, NTR 
IF  (  ISEP  .EQ.  1  )  WRITE  (6,110)  XCTR, XTR, NTR 
ICODE   =  1 


C 
C 


f**  CALCULATE  INTERMITTENCY  DISTRIBUTION  *** 

RXNTR  =  XTR  *  UETR  *  RL 

GGFT   =  RL**2/PGAMTR/RXNTR**1.  34*UETR**3 

DO  30  I=NTR1,NXT 

ALFAS(I)  =  0. 0168 

GMTRS(I)=  1.  0 
30  CONTINUE 

ALFAS(NTR)  =  0. 0168 

UEINTG  =  0.0 

Ul  =  0.5/UETR 

XI  =  XTR 

DO  40  I=NTR1,NXT 

U2  =  0. 5/UE(I) 

X2  =  X(I) 

UEINTG  =  UEINTG+(U1+U2)*(X2-X1) 
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Ul  =  U2 
XI   =  X2 

GG  =  GGFT*UEINTG*(X(I)-XTR) 
IF(GG    . GT.     10. 0)    GOTO   50 
GMTRS(I)    =    1. O-EXP(-GG) 
40   CONTINUE 


C 

c 


60 

C 
C 

C 

C 


10 

20 


30 


C 
C 
C 
C 
C 
C 
C 
C 
C 


f**  RESET  FINITE  DIFFERENCE   CALCULATIONS  *** 
50  DO   60   J=1,NPT 

F(J,2)   =  F(J,1) 

U(J,2)   =  U(J,1) 

V(J,2)   =  V(J,1) 

B(J,2)   =  B(J,1) 

W(J,2)   =  W(J,1) 

CONTINUE 

RETURN 

END 


SUBROUTINE   EDGCHK(NP,    ETA,    F,    U,    V) 

DIMENSION  ETA(lOl),    F(101),    U(101),    V(101) 

JS  =  NP    -    3 

NPM1      =  NP    -    1 

DO    10   J=JS,    NPM1 

JJ  =   J 

IF(U(J).GE.U(NP)    .OR.    V(J).LT.  0.  0)    GOTO   20 

CONTINUE 

RETURN 

JS  =  JJ   -    1 

IF(JS.GT.  (NP-2))    JS  =  NP-2 

CALL  AMEAN(JS,    NP,    ETA,    U,    1) 

CALL  AMEAN(JS,    NP,    ETA,    F,    1) 

DETAP     =  ETA(JS)    -ETA(JS-l) 

VJP  =   (U(JS)-U(JS-l)) /DETAP 

DO   30   J=JS,NPM1 

DETAM     =  ETA(J+1)-ETA(J) 

VJM  =   (U(J+1)-U(J))/DETAM 

V(J)        =   (VJM*DETAP  +  VJP  *  DETAM)  /(DETAP+DETAM) 

VJP  =  VJM 

DETAP      =  DETAM 

CONTINUE 

V(NP)      =   -V(NP-l)   +2.0  *  VJP 


RETURN 


*-  .*-  -•-  **-  -'-  -'-  -'-  »'-  J-  -»-  J-  J-  -'-  JL  »'-  -*-  _'-  -J-  -'-  -'-  J1-  J-  «'-  «*-  J-  J"-  .'-  «'*  «.»-  J-  A  .'-  JL  J-  JU  «»-  JL  ..»,  -»,  -T.  JL,  Jm  .!.  J*  JU  <Jm  -'-  -J-  - 

'%   t\    t\    -  »    /  .    .  <■    *\   *>   n   **  t\    '  v    st    /t    .  t    /t   ft   #*    r*    /»   /V/v   *»    /\    #v   *\    0\    **   **    #\   /W5/W*   **   **   *V  /»   **   *\   #»   #*   m7v  *\   *\   *\   *»    I 


i  '  '  »V  **  «t  rfv  # 


NOTES: 


(FOR  CHANGING  FROM  THE   ORIGINAL  PROGRAM) 


4t 


1.  'EDDY'    HAS   BEEN  MODIFIED  BY  ADDING    'FINT'. 

2.  SUBROUTINE    'EDGCHK'    HAS   BEEN  ADDED. 

3.  GROWTH  LIMIT  HAS   BEEN  ADDED  FOR  2   IN    'MAIN2'. 


* 


END 
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SUBROUTINE  SMPSON 


C 

c 
c 

c 
c 
c 
c 
c 
c 

C( 

c 
c 
c 
c 
c 
c 
c 
c 

C( 

c 


30 


35 
38 


c 
c 


+ 


COMMON  /BLCO/ 
COMMON  /BLC1/ 
COMMON  /BLC7/ 
COMMON /EDDY 1/ 

COMMON  /GTY  / 
COMMON  /GRD  / 
DIMENSION  CRD( 
DATA  RTD/O. 00, 
f- 
DATA  CRD/5. 00, 


+ 


NX,NXT,NP,NPT,NTR,IT,INVRS,NS,IP 
F(101,2),U(101,2),V(101,2),W(101,2),B(101,2) 
C(100,100),D(100),DB(100),DBP(100) ,UE0(100),GI 
RL , RX , SQRX , RXNTR , GMTR , GMTRS (100) 

, ALFAS( 100) ,FFS( 100) ,RTS( 100) , IEDY,NXSPT 
X(101),UE(100),P1(100),P2(100),CEL,CELH 
ETA( 101) ,DETA( 101) , A( 101) 
12),RTD(12) 

0.  05,0.  12,0.  20,0.  30,0.  40,0.  50,0.  60,0.  70, 
0.  80,0.  90,1.  00/ 

4.  75,4.  35,3.  80,3.  25,2.  85,2.58,2.  37,2.  25, 
2.  15,2.  06,2.  00/ 


STEP  1  CALCULATE  (DU/DX)/(DU/DY) 
IF(NX.  LT.  NXSPT)  GOTO  10 

IN  THE  SEPARATED  REGION,  ALFA  SET  TO  BE  CONSTANT 
ALFAS(NX)=  ALFASP 
RETURN 


10  CONTINUE 

IF(V(1,2).GT.  0.  0)  GOTO  20 

SEPARATION  OCCURS.  ALFA  SET  TO  BE  THE  PREVIOUS  ITERATED  VALUE 
ALFASP  =  ALFAS(NX) 
NXSPT  =  NX 
RETURN 


MODIFY  OUTER  EDDY  BASED  ON  SIMPSON  SUGGESTION 
TM     =0.0 
JM     =  1 
DO  30  J=2,NP 

TS     =  (B(J,2)-1.0)*  V(J,2) 
IF(TS.  LT.TM)  GOTO  30 
TM     =  TS 
JM     =  J 
CONTINUE 
VNXM  =  0.5*(V(JM,2)+V(JM,1)) 

IF  (NX  .LE.  NS)  GOTO  35 

DUDX  =  (U(JM,2)-U(JM,1))  /  (X(NX) -X(NX-l) ) 
GO  TO  38 

DUDX   =  CEL*(U(JM,2)-U(JM,1))+P2(NX)*U(JM,2)+0.5*ETA(JM)* 
+  VNXM*(P2(NX)-1.  0) 

RU     =  RL 

IF(NX. LE.NS)RU  =  RL  *  UEO(NX)  *  X(NX) 
RL2    =  SQRT(RU) 
RR     =  DUDX/VNXM/RL2 

STEP  2  :  CALCULATE  (UU  -  VV)/UV 
VNXM  =  0.5*(V(1,2)+V(1,1)) 
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c 
c 
c 
c 
c 
c 

60 

c 

c 

70 


RT    =  VNXM/TM 
PRINT'(3X,2I5,3F10.3)' ,NX, JM, VNXM,TM,RT 

IF  (RT  .  LT.  0.  0)  RT  =  0.  0 

IF(RT.  GT.  1.  0)  GOTO  60 

CR    =6.0  /(1.0  +   2.0  *  RT*(2.  0  -RT)) 

CR    =2.0 

DO  40  1=2,12 

IF(RT.  LT.  RTD(I))  GOTO  50 
40  CONTINUE 

GOTO  70 
50  CR    =  CF.i;i-i)+(CRD(I)-CRD(I-l))*(RT-RTD(I-l))/(RTD(I)-RTD(I' 

GOTO  70 

CR    =  ( 1.  0  +  RT  )  /RT 

STEP  3  :  CALCULATE  FF 
FR   =  CR  *  rr 
IF(FR  .GT.  0.  35)  FR  =  0.  35 
IF  (FR  .LT.  -0.  8)  FR  =  -0.  8 
FFS(NX)=  (FFS(NX)  +   (1.0  -FR))/  2.0 
RTS(NX)=  RT 

ALFAS(NX)=  0. 0168/FFS(NX)**2.  5 
RETURN 


C( 
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20 
30 


40 


END 

SUBROUTINE  XSPACE(NI , NRITE , XII ,XLLT, RAD, NL1,NR1) 

DIMENSION  XII(200),T(200) 

DATA  PI/3. 14159265359879/ 

RAD  =  PI 

NLEFT=NI-1-NRITE 

NR4=NRITE/2 

IF((NRITE/2*2)  .  NE.  NRITE)  NR4=(NRITE+l)/2 

NL1=NR4+1 

NL2=NR4+NLEFT 

NR1=NL2+1 

NR2=NI 

PI 2=0.  5*PI 

RAD2=(PI-RAD)/2. 0+PI2 

RAD3=RAD2+RAD 

SRT  =RAD2/FL0AT(NR4) 

SRT2=SRT 

IF(( NRITE/ 2*2)  .  NE.  NRITE)  SRT2=RAD2/FL0AT(NR4-1) 

SLT  =  RAD/FLOAT(NLEFT) 

DO  10  1=1, NR4 

XII(I)=0.5*(1.0+COS(FLOAT(I-1)*SRT)) 

DO  20  I=NL1,NL2 

XII(I)=0.  5+XLLT*COS(FLOAT(I-NLl)*SLT+RAD2) 

DO  30  I=NR1,NR2 

XII(I)=0.5*(1.0+COS(FLOAT(I-NR1)*SRT2+RAD3)) 

NA=(NI+l)/2 

IF((NI/2*2)  .EQ.  NI)  NA=NI/2+l 

FN1=FL0AT(NA-1) 

FN2=FN1 

IF((NI/2*2)  .EQ.  NI)  FN2=FLOAT(NA-2) 

DO  40  1=1, NA 

T(I)=FL0AT(NA-I)/FN1 
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45 


50 


C 
C 


CALL  AMEAN(1,NA,T,XII,1) 
XDIF  =  XII(l)  -  XII(2) 
IF(XDIF  . LT.  0. 004)  THEN 
DO  45  1=2,5 

XI  1(1)  =  XII(I-1)-XDIF*3.  0 
CONTINUE 

CALL  AMEAN(2,NA,T,XII,10) 
END  IF 

DO   50    I=NA,NI 
XII(I)   =  XIKNI-I+1) 
RETURN 
END 

SUBROUTINE  TRGRID  (Nl  ,  XO  . 
THIS  SUB.  IS  TO  REGRID  SPACING 


C 

c 
c 
c 
c 


60 


Y0,NI,NRITE,XLLT,N10,RAD,ID,NXSS) 
NEAR  TRAILING -EDGE 


62 
65 

68 


DIMENSION  X0(200) ,YO(200) ,XI(200) ,YI(200) ,D1(200) ,D2(200) ,D3(200) 
+  XOOf 200) ,YOO(200) ,XII(200) ,YII(200) ,WX(200) ,WY(200) , 

+  WXI  200),WYI(200),T(200) 


N20    =  N10 
IF((Nl/2*2)  .EC 
IF((NI-(NJ/2)''-2) 

IFUNI-(NI/2)*2) 
N2I    =  Nil 
IF((NI/2*2)  .EQ. 


Nl)  N20  =  N10+1 

. NE.  0)  N1I=  (NI-l)/2+l 

.EQ.  0)  NlI=NI/2 

Nl)  N2I  =  N1I+1 


CALL  XSPACE(NI,NRITE,XI,XLLT,RAD,NL1,NR1) 

PRINT  * , ' NRITE= ' , NRITE , '  XLLT= ' , XLLT 
WRITE  (6,  290) 

WRITE  (6,  300)  (XO(I)  ,1=1, Nl) 
WRITE  (6,  298) 

WRITE  (6,  300)  (YO(I)  ,1=1, Nl) 
IF(ID  .EQ.  2)  THEN 
DO  60  I=NL1,NXSS 
YI(I)=Y0(I) 
XI(I)=X0(I) 
NXST=NXSS+8 

XMl=(XI(NXST)-XI(NXSS))/8.  0 
XM2=(XI(NR1)-XI(NXSS))/(NR1-NXSS) 
DO  62  I=NXSS,NR1-1 
XI(I)=(XO(I)+XI(I))/2.  0 
DO  65  I=NXSS,NXST 
T(I)=FL0AT(I-NXSS)*XM1+XI(NXSS) 
CALL  AMEAN(NXSS,NXST,T,XI,8) 
DO  68  I=NXSS,NR1 
T(I)=FLOAT(I-NXSS)*XM2+XI(NXSS) 
CALL  AMEAN(NXSS,NR1,T,XI,12) 
N10=NL1 
N1I=N10 
N20=N1+1-NXSS 
N2I=NI-NXSS+1 
ELSE 

NXSS=N2I 
END  IF 
FOR  LOWER  SURFACE 
DO  5  I  =  1  ,  N10 
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9 
C 
C 


10 


20 


II  =  N10  -  I  +  1 

VX(II)  =  XO(I) 

WY(II)  =  YO(I) 

CONTINUE 

DO  7  I  =  1  ,  Nil 

II  =  Nil  -  I  +  1 

WXI(II)  =  XI(I) 

CONTINUE 

CALL  DIFFSCNIOjWXjWYjDI.DZ^jO) 

CALL  INTRP3(N10,WX,WY,D1,D2,D3,N1I,WXI,WYI) 

DO  9  I  =  1  ,  Nil 

II  =  Nil  -  I  +  1 

YI(II)  =  WYI(I) 

CONTINUE 

FOR  UPPER  SURFACE 
DO  10  I  =  1  ,  N20 
II  =  Nl  -  N20  +  I 
XOO(I)  =  XO(II) 
YOO(I)  =  YO(II) 
CONTINUE 

DO  20  I  =  1  ,  N2I 
II  =  NI  -  N2I  +  I 
WXI(I)  =  XI(II) 
CONTINUE 

CALL  DIFF3(N2O,XOO3YOO,Dl,D2,D3,0) 
CALL  INTRP3(N20,XOO,YOO,Dl,D2,D3,N2I,WXI,WYI) 
DO  25  I  =  1  ,  N2I 
II  =  N2I  -  I  -■■  1 
XII(II)  =  WXI(I) 


YII(II)  = 
CONTINUE 


VYI(I) 


25 
C 

C   COMBINE  TOO  SURFACES  INTO  ONE  CIRCLE 

C 

C     NN  =  Nl  -  N10  -  N20 

C 

C 

C 

C 

C 

C30 


NN 


N10 

1  , 
+  I 
+  I 

XO(III) 
YO(III) 


40 


NN  =  Nl  - 

DO  30  I  = 

II  =  Nil 

111=  N10 

XI(II)  = 

YI(II)  = 

CONTINUE 

DO  40  I  =  1  ,  N2I 

11  =  NXSS  +1-1 

12  =  N2I  -  I  +  1 
XI(I1)  =  XII(I2) 
YI(I1)  =  YII(I2) 
CONTINUE 

XI(1)  =  XO(l) 
XI(I1)=  XO(Nl) 
YI(1)  =  YO(l) 
YI(I1)=  YO(Nl) 


Nl  =  II 
DO  50  I 
XO(I)  = 


=  1  , 
XI(I) 


Nl 
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50 

C 

C 

C 

C 

c 
c 


YOfI)  =  YI(I) 
CONTINUE 


C  - 
100 

200 
290 
295 
298 
300 
C  - 


c 
c 
c 


10 

c 


c 
c 
c 
c 
c 
c 
c 


20 


WRITE  (6 
WRITE  (6 
WRITE  (6 
WRITE  (6 

RETURN 


295) 

300)  (XO(I)  ,  1=1, Nl) 

298) 

300)  (YO(I)  ,  1=1, Nl) 


FORMAT(7I5) 

FORMAT(6F10.  0) 

FORMATC / , '   ORIGINAL  COORDINATES ' , / , '   X/C  * ) 

FORMAT( / , '   INTERPLATED  COORDINATES ' , / , '   X/C  '  ) 

FORMATC      Y/C') 

FORM AT (6F 10.  6) 


END 

SUBROUTINE  STAGR( N , STAG , XO , YO , XSTGR , YSTGR ) 

DIMENSION  XO( 100) ,YO( 100) ,XSTGR( 100) ,YSTGR( 100) ,DS( 100) 

XOTE  =  0.5  *  (XO(l)+XO(N)) 

YOTE  =  0.5  *  (YO(l)+YO(N)) 

DS(1)  =  SQRT((XO(l)-XOTE)**2  +  ( YO( 1) -YOTE)**2) 

DSM  =  DS(1) 

DO  10  I  =  2  ,  N 

DS(I)  =  SQRT((XO(I)-XOTE)**2  +  (YO(  I)  -YOTE)--* 2) 

IF  (DS(I)  .  LT.  DSM)  GOTO  10 

IM  =  I 

DSM  =  DS(I) 

CONTINUE 

YYY  =  YOTE-YO(IM) 

XXX  =  XOTE-XO(IM) 

IF  (YYY  .EQ.  0.0  .AND.  XXX  .  EQ.  0.0)  THEN 

ANG  =  0.0 

ELSE 

ANG  =  ATAN2(Y".r,XXX) 

END  IF 

ANG  =  ANG  +  STAG 

COS AN  =  COS (ANG) 

SINAN  =  SIN(ANG) 

DO  20  I  =  1  ,  N 

YY  =  YO(I)-YO(IM) 

XX  =  XO(I)-XO(IM) 

IF  (YY  .EQ.  0.0  .AND.  XX  .  EQ.  0.0)  THEN 

ANGCO  =  0.0 

ELSE 

ANGCO  =  ATAN2(YY,XX) 

END  IF 

XSTGR(I)=  XO(I)*COSAN  +  YO(I)*SINAN 

YSTGR(I)=  YO(I)*COSAN  -  XO( I)*SINAN 

CONTINUE 
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RETURN  INT30070 

END  INT30080 
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APPENDIX  B.     C4  CASCADE 

A.     EXPERIMENTAL  RESULTS 

The  experimental  results  of  the  C4  cascade  were  obtained  directly  from  professor 
G.J.  Walker,  University  of  Tasmania,  Tasmania,  Australia,  who  performed  these  exper- 
iments. 

The  results  of  the  boundary  layer  measurements  of  the  C4  cascade  are  given  below 
at  four  inlet  angles:  34.1°,  36.3°,  45.6°,  and  47.7°.  The  Reynold  numbers,  based  on  the 
chord  and  the  upstream  velocity,  are  200000,  191000,  173000  and  171000  respectively. 
The  results  given  in  the  following  tables  include  the  displacement  thickness  (<5'),  the 
shape  factor  (H)  and  the  local  free  stream  velocity  (UE). 

Table  1.     EXPERIMENTAL  RESULTS  AT  INLET  ANGLE  OF  34.1° 


X  c 

6*  [10"3FT] 

H 

UE  [FT  SEC] 

0.4 

4.9 

2.4S 

168.37 

0.5 

6.28 

2.61 

167.35 

0.6 

8.79 

3.24 

158.31 

0.7 

10.83 

3.63 

149.27 

O.S 

16.63 

3.79 

147.13 

0.9 

16.19 

1.S9 

143.79 

Table  2.     EXPERIMENTAL  RESULTS  AT  INLET  ANGLE  OF  36.3C 


X  c 

(M10~4FT] 

H 

UE  [FT  SEC] 

0.4 

5.43 

2.55 

161.63 

0.5 

7.09 

2.70 

157.59 

0.6 

10.3 

3.34 

14S.7S 

0.7 

12.63 

3.78 

139.S7 

0.8 

14.84 

2.78 

135.01 

0.9 

16.43 

1.76 

133.23 
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Table  3.     EXPERIMENTAL  RESULTS  AT  INLET  ANGLE  OF  45.6C 


X  c 

6*  [10~3  FT] 

H 

UE  [FT  SEC] 

0.4 

8.08 

2.58 

137. 8S 

0.5 

9.83 

2.41 

133.70 

0.6 

12.35 

2.33 

122.18 

0.7 

12.9S 

1.97 

114.93 

0.8 

19.44 

1.90 

111.77 

0.9 

27.69 

1.92 

109.26 

Table  4.     EXPERIMENTAL  RESULTS  AT  INLET  ANGLE  OF  47.7° 


X  c 

c5*[10-4FT] 

H 

L'E  [FT  SEC] 

0.4 

8.87 

2.24 

130.64 

0.5 

10.27 

2.19 

124.76 

0.6 

14.31 

2.08 

116.86 

0.7 

16.45 

1.87 

106.72 

0.8 

24.16 

1.82 

103.75 

0.9 

36.00 

2.01 

102.18 
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The  results  of  the  measurements  of  the  velocity  profiles  in  the  boundary  layer  at  two 
inlet  aneles.  34.1°  and  36. 3C  at  50%  chord  are  given  below. 


Table  5.     VELOCITY  PROFILES  AT  50%  CHORD. 


V 

/?  =  36.3° 

/?  =  34.1° 

0.0 

0.0 

0.0 

2.3 

0.172 

0.20S 

3.7 

0.270 

0.327 

6.2 

0.469 

0.534 

8.6 

0.666 

0.728 

1 1.0 

0.794 

0.867 

13.4 

0.S91 

0.933 

18.3 

0.982 

0.985 

23.2 

1 .000 

1.000 

B.     C4  CASCADE  COORDINATES 


DIMENSION  X(0:  100),XU(0:  100),XL(0:  100),YU(0:  100),YL(0:  100) 
DATA  Al,A2,A3,A4/0. 15492,0. 06563,0. 2528,0. 2811/ 
DATA  Bl,B2,B3,B4/0. 03866,0. 07871,0. 1467,0. 03448/ 
PI  =  AC0S(-1.  0) 
C     READ  (5,800)  NMAX 
800  FORMAT  (15) 
NMAX=33 
C     READ  (5,810)  (X( I) , 1=0 ,NMAX) 
810  FORMAT  (6F10. 6) 
DO  50  1=0, NMAX 

X(I)  =  (1. 0-C0S(PI*I/NMAX))/2. 
50  CONTINUE 

DO  100  1=0, NMAX 

SRT  =  SQRT((0.5/SIN(PI/12))**2-(0.5-X(I))**2) 

YC  =  -0.5/TANCPI/12)  +  SRT 

DY  =  ATAN((0.5-X(I))/SRT) 

IF  (X(I).LT.  0.  3)  THEN 

YT  =  A1---SQRT(X(I))  -  A2*X(I)  - 
ELSE 

YT  =  Bl  +  B2»X(I)  -  B3*X(I)**2 
END  IF 


A3*X(I)**2  +  A4*X(I)**3 
+  B4*X( I)** 3 


C4 
C4 
C4 
C4 
C4 
C4 
C4 
C4 
C4 
C4 
C4 
C4 
C4 
C4 
C4 
C4 
C4 
C4 
C4 
C4 
C4 


00010 
00020 
00030 
00040 
00050 
00060 
00070 
00080 
00090 
00100 
00110 
00120 
00130 
00140 
00150 
00160 
00170 
00180 
00190 
00200 
00210 


YU(I)  =  YC  +  COS(DY)*YT 

YL(I)  =  YC  -  COS(DY)*YT 

XU(I)  =  X(I)  -  SIN(DY)*YT 

XL(I)  =  X(I)  +  SIN(DY)*YT 
100  CONTINUE 
C     WRITE  (6,900)  ( I ,X( I) ,XU( I) , YU( I) ,XL( I) ,YL( I) , I=0,NMAX) 
C  900  FORMAT  ( 15 ,4X ,F10. 6,4X,2F10. 6 ,4X,2F10.  6) 

WRITE  (1,910)  (XL(I),I=NMAX,0,-1),(XU(I),I=1,NMAX) 

WRITE  (1,910)  (YL(I),I=NMAX,0,-1),(YU(I),I=1,NMAX) 
910  FORMAT  (6F10.6) 

STOP 

END 


C4  00220 
C4  00230 
C4  00240 
C4  00250 
C4  00260 
C4  00270 
C4  00280 
C4  OO290 
C4  OO3O0 
C4  00310 
C4  00320 
C4  OO330 
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